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This  study  explores  the  optimization  method  of  simulated  annealing  for  use  in  an 
analytic  tool  for  counter-narcotics  analysis.  A  model  is  developed,  employing  RAND 
Corporation's  Simulation  of  Adaptive  Response  model  as  an  objective  function  evaluator, 
which  optimizes  interdiction  asset  locations  relative  to  a  sample  smuggling  network  in  the 
Caribbean  region.  In  addition  to  asset  location  optimization,  the  response  of  the  model 
to  changing  numbers  of  assets  is  also  tested.  Results  indicate  that  this  methodology  has 
potential  for  use  in  the  counter-narcotics  program,  and  perhaps  other  network  interdiction 
applications.  Further  research  and  testing  are  recommended.  Military  drug  interdiction 
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I.    INTRODUCTION 

A.   BACKGROUND 

The  Defense  Authorization  Act  of  1989  designated  the  Department  of  Defense 

(DoD)  "as  the  single  lead  agency  of  the  Federal  Government  for  the  detection  and 

monitoring  of  aerial  and  maritime  transit  of  illegal  drugs  into  the  United  States."  The  Act 

also  tasks  DoD  with  expanded  support  to  civilian  Law  Enforcement  Agencies  (LEA's) 

involved  in  drug  interdiction.   According  to  Secretary  of  Defense  Dick  Cheney, 

The  Department  of  Defense  is  an  enthusiastic  participant  in  the  nation's  drug 
control  effort.  We  have  significant  resources  at  our  disposal.  We  can  make  a 
substantial  contribution  to  our  national  effort  if  we  use  our  assets  intelligently  and 
efficiently.  [Ref.  1:  p.  41 

To  these  ends,  the  President's  National  Drug  Control  Strategy  for  1991  estimated  the 

DoD  budget  for  drug  interdiction  at  nearly  $1.1  billion,  more  than  double  the  actual 

expenditures  for  1989.    Obviously,  the  potential  contributions  by  DoD  are  being  taken 

very  seriously. 

In  1988,  the  RAND  Corporation  developed  the  Simulation  of  Adaptive  Response 

(SOAR)  model  to  analyze  the  effect  of  military  involvement  in  drug  interdiction.  SOAR 

is  a  stochastic  network  simulation  of  drug  smuggling  involving  theoretical  routes  for  drug 

shipments  and  the  response  of  smugglers  to  being  interdicted.  Expanded  involvement  of 

the  military  in  drug  interdiction  was  evaluated  by  examining  the  effect  that  increased 

interdiction  probability  on  the  theoretical  routes  might  have  on  the  import  price  of  illegal 


drugs.  The  conclusions  of  the  authors  of  the  RAND  study  (Gordon  Crawford,  Peter 
Reuter,  et  al.)  were  not  heartening.  Although  the  limitations  of  the  SOAR  model  caused 
it  to  be  biased  toward  successful  interdiction,  it  was  found  that  "interdiction  must  be  very 
stringent  indeed  to  greatly  affect  U.S.  drug  consumption."  [Ref.  2:  p.  73] 
However,  the  model  is  simplistic  in  its  approach  to  smuggling  routes  and  the  presence  of 
interdiction  forces.  Routes  are  not  associated  with  actual  geographic  areas  due  to  a  lack 
of  available  data,  while  interdiction  assets  are  represented  solely  by  a  probability  of 
interdiction  on  the  generic  routes.  This  simplicity  leaves  room  for  expansion  of  the  model 
and  further  analysis. 

B.       PROBLEM  ANALYSIS 

Government-sponsored  studies  prior  to  that  undertaken  by  RAND  did  not 
concentrate  on  the  adaptation  of  smugglers,  nor  explicitly  examine  the  potential  for 
military  contribution  to  the  drug  interdiction  program  [Ref.  3:  p.  4].  As  part 
of  a  comprehensive  study,  SOAR  was  designed  to  aid  in  a  policy  analysis  of  the  overall 
interdiction  program. 

As  a  policy  tool,  SOAR  did  not  require  a  geographic  structure  or  specific 
representation  of  interdiction  assets.  For  this  reason,  SOAR  is  not  useful  in  operational 
planning  or  analysis.  RAND's  conclusions,  at  a  policy  level,  indicate  that  military 
involvement  in  drug  interdiction  may  have  limited  effectiveness  [Ref.  3:  p.  130]. 
However,  the  Department  of  Defense  has  a  major  role  in  the  program,  and  would  benefit 
from  the  availability  of  an  analytic  tool  which  could  be  used  to  evaluate  the  effects  of 


operational  planning  and  to  explore,  with  a  consistent  measure,  the  effectiveness  of 
different  force  and  capability  alternatives. 

Interdiction  problems  are  frequently  modeled  as  flow  networks  for  which 
optimization  methods  are  used  to  determine  which  arcs  or  nodes  should  be  cut  or 
removed  to  provide  an  optimum  result.  In  the  case  of  a  smuggling  network,  transit  routes 
may  be  represented  by  arcs  which  are  not  well  defined,  or  are  completely  unknown  to  the 
interdictors,  making  it  difficult  to  entirely  sever  the  arc.  Even  if  it  were  possible,  the 
smugglers  themselves  will  eventually  realize  that  the  arc  has  been  cut  and  will  move  to 
another  route,  adapting  their  methods  to  the  interdiction.  The  SOAR  model  simulates  this 
type  of  response  on  the  part  of  drug  smugglers,  making  it  a  fairly  obvious  choice  to 
measure  the  effectiveness  of  interdiction  on  a  smuggling  network. 

Because  routes  used  by  smugglers  may  not  be  clearly  defined,  it  is  unrealistic  to 
merely  place  assets  on  arcs  representing  segments  of  the  routes.  By  developing  a 
geographic  structure  for  the  routes,  assets  can  be  located  by  latitude  and  longitude,  and 
the  effects  of  assets  on  the  routes  (i.e.,  the  probability  of  interdiction  on  routes)  can  be 
measured  as  a  function  of  their  range  from  the  routes.  This  also  would  allow  assets  to 
affect  multiple  routes,  for  example,  routes  passing  through  a  geographically  restricted 
area,  or  chokepoint.  Finding  the  optimum  locations  for  assets  relative  to  a  smuggling 
network  with  stochastic  flows  will  be  difficult.  As  the  locations  of  assets  change,  flows 
through  the  routes  will  also  change,  resulting  in  local  and  global  optima  which  vary  with 
the  locations  of  assets.  There  are  also  likely  to  be  many  such  optima,  a  result  of  different 
assets  contributing  the  same  probability  of  interdiction  to  routes  or  being  located  in 


positions  previously  tested  by  other  assets.  In  addition,  the  number  of  possible  locations 
for  assets  is  infinite,  requiring  some  control  over  the  locations  tested  if  discrete 
optimization  is  to  be  used.  This  control  must  include  geographically  feasible  areas,  but 
should  also  reflect  operational  constraints,  such  as  a  requirement  for  certain  assets  to 
operate  within  a  set  range  from  each  other  or  a  base,  and  political  constraints,  such  as 
foreign  airspace  or  waters  restrictions. 

An  optimization  model  with  the  properties  described  in  this  analysis,  which  retains 
enough  realism  to  be  considered  valid,  would  be  of  great  benefit  at  the  operational  level 
of  the  military  drug  interdiction  program. 

C.   METHODOLOGY 

We  will  propose  a  solution  to  the  need  for  an  analytical  tool  which  may  be  used  for 
analysis  and  planning  at  the  operational  level  of  the  military  drug  interdiction  program. 
An  exploration  of  the  use  of  simulated  annealing  to  optimize  the  location  of  counter- 
narcotics  assets  with  respect  to  a  drug  smuggling  network  is  presented.  The  objective 
function  is  an  outcome  of  the  modified  SOAR  model,  in  which  the  routes  used  by  the 
smugglers  are  modeled  with  a  geographical  structure,  and  the  physical  location  of 
interdiction  assets  (eg.  ships,  aircraft,  aerostat  balloons,  observation  posts,  etc.)  are  used 
to  determine  the  probability  of  interdiction  in  an  area.  From  the  SOAR  model  we  extract 
the  total  cost  to  smugglers,  a  measure  of  effectiveness  for  one  major  goal  of  the  overall 
drug  interdiction  program,  reduction  in  consumption  of  illegal  drugs.  Interdiction  raises 
the  risks  involved  with  drug  smuggling,  thereby  increasing  the  costs  associated  with 
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personnel  pay,  equipment  and  drug  replacement,  etc.  Eventually,  this  results  in  higher 
consumer  prices  for  drugs.  Economically,  higher  prices  should  cause  a  decrease  in 
demand  for  drugs  in  the  long  run,  that  is,  lower  consumption  by  Americans.  Increased 
crime  among  addicts  to  compensate  for  higher  prices  is  not  addressed. 

Use  of  SOAR  as  the  objective  function  evaluator  in  a  model  which  exhibits 
numerous  local  and  global  optimal  asset  locations  eliminates  the  potential  use  of  gradient- 
optimization  methods.  It  is  for  this  reason  that  we  chose  to  apply  simulated  annealing 
as  the  optimization  method  to  maximize  the  total  cost  to  the  smugglers. 

D.       SCOPE 

This  thesis  briefly  discusses  the  flow  of  cocaine  from  South  and  Central  America 
into  the  southern  United  States,  and  the  role  the  Department  of  Defense  plays  in 
interdiction  efforts  within  this  geographical  area. 

The  model  presented  is  capable  of  evaluating  the  use  of  assets  in  the  Atlantic, 
Pacific  and  Southwest  Border  regions;  however,  for  the  numerical  results  presented  in  this 
thesis,  only  routes  and  assets  in  the  Caribbean  and  Gulf  of  Mexico  are  examined.  Much 
of  the  information  that  is  available  about  known  smuggling  routes  and  the  assets  involved 
in  counter-narcotics  operations  is  highly  classified.  In  the  interest  of  maintaining  an 
unclassified  example,  some  data  used  in  the  model  is  taken  from  the  original  SOAR 
model,  while  information  on  smuggling  routes  and  counter-narcotics  assets  is  intended  to 
be  representative  only,  and  should  not  be  construed  as  data  representing  actual 
circumstances. 


II.     THREAT  AND  INTERDICTION  EFFORTS  OVERVIEW 

A.       CURRENT  THREAT  OF  ILLEGAL  DRUG  TRAFFICKING 

A  number  of  indicators  show  that  progress  has  been  made  in  the  war  against  drugs. 
Overall  drug  use,  drug-related  medical  emergencies,  and  student  attitudes  toward  drug  use 
have  all  declined  past  the  goals  set  in  the  1989  National  Drug  Control  Strategy,  shown 
in  Table  1.  [Ref.  4:  p.  5] 


Table  1     NATIONAL  DRUG  CONTROL  STRATEGY  NATIONAL  OBJECTIVES:     Statistical 
significance  ranges  from  1  to  30%. 


GOAL 

ACTUAL 

CURRENT  OVERALL  DRUG  USE 

-  10% 

-  11% 

CURRENT  ADOLESCENT  DRUG  USE 

-  10% 

-  13%' 

OCCASIONAL  COCAINE  USE 

-  10% 

-  29% 

FREQUENT  COCAINE  USE 

50% 

reduction  in 

rate  of  increase 

-  23% 
actual  decline 

CURRENT  ADOLESCENT  COCAINE  USE 

-  20% 

-49% 

DRUG-RELATED  MEDICAL  EMERGENCIES 

-  10% 

-  18% 

STUDENT  ATTITUDES  TOWARD  DRUG  USE 

-  10% 

-  28% 

Cocaine-specific  indicators  showed  decreased  purity  and  significantly  higher  prices 
in  many  areas  of  the  country  for  the  past  year.  Also,  since  1985,  "...  reported  current 
cocaine  use  (use  in  the  past  month)  decreased  some  72  percent." 
[Ref.  5:  p.  1] 


The  effect  on  smugglers  themselves  may  be  deduced  from  the  surrender  of  a 
number  of  Colombian  drug  lords  and  the  level  of  effort  smugglers  put  forth  to  get  their 
product  into  the  United  States.  Examples  of  this  effort  include  the  tunnel  between 
Mexico  and  Arizona  discovered  in  May  1990  [Ref.  6:  p.  Al]  and  the 
occasionally  high  tech  hiding  places  used  in  vehicles,  boats,  and  aircraft. 

On  the  other  hand,  many  problems  still  exist  which  illuminate  the  threat  to  national 
security  declared  by  President  Reagan  in  1986.  Drug-related  violence  and  murders  remain 
at  high  levels  within  the  United  States,  as  does  the  amount  of  cocaine  being  produced 
throughout  South  America.  Table  2  displays  production  levels  for  the  three  major 
producing  countries  from  1988  to  1990,  the  last  year  for  which  data  was  available. 


Table  2    COCAINE  PRODUCTION  ESTIMATE:    Potential  Cocaine  HC1  Production  by  Country 
(tons).    Sources:   NNICC  Reports,  1989  and  1990. 


COUNTRY 

YEAR 

1988 

1989 

1990 

PERU 

242 

440-468 

440-473 

BOLIVIA 

112 

264-458 

259-435 

COLOMBIA 

43 

77 

72 

TOTAL 

397 

781-996 

771-980 

Smugglers'  high  tech  methods  display  their  increasing  adaptability  to  interdiction 
efforts,  as  well  as  the  effects  of  interdiction.  Some  smugglers  are  showing  that  they  are 
willing  (or  crazy  enough)  to  take  on  interdiction  forces,  as  evidenced  by  incidents  along 
the  border  with  Mexico  where  smugglers  have  exchanged  small  arms  fire  with  U.S. 


military  forces  [Ref.  7:  p.  1],  and  the  attempt  in  April  1990  by  smugglers 
linked  to  the  Medellin  cartel  to  buy  Stinger  missiles  [Ref.  8:  p.  3]. 

The  evolution  of  narco-terrorism  may  be  proof  enough  that  we  are  indeed  facing 
a  war  which  must  be  fought  by  our  armed  forces  outside  the  United  States.  Many 
revolutionary/terrorist  organizations,  such  as  M-19  and  FARC  in  Colombia  and  Shining 
Path  in  Peru,  are  intimately  linked  to  the  drug  trade,  forming  the  basis  for  this  narco- 
terrorism.  These  groups  have  shown  decided  anti-American  leanings,  and  have  used 
violence  and  threats  against  high  level  U.S.  officials.  One-time  M-19  leader  Ivan  Mariano 
Ospina  said  in  December  1984,  "May  these  threats  be  carried  out  and  may  they  be  carried 
out  in  the  entire  world  against  the  rapacious  imperialism  that  lives  at  the  cost  of  misery 
of     exploited      people...."      [Ref.  9:  pp.  36,37]  Even      state      leaders      have 

been  involved,  most  recently  illuminated  by  General  Noriega  in  Panama,  but  also 
including  Fidel  Castro  who  has  acknowledged  the  use  of  cocaine  as  a  weapon,  "We  are 
going  to  make  the  people  up  there  (the  United  States)  white,  white  with  cocaine  [Ref.  9: 
P-  32]." 

Clearly,  there  are  some  signs  of  success  in  the  drive  to  reduce  illegal  drug 
trafficking  and  abuse.  However,  it  is  also  plain  that  much  more  effort  is  required  in  all 
phases  of  the  war  on  drugs  before  the  situation  can  be  considered  under  control,  let  alone 
resolved.  Involvement  of  the  Department  of  Defense  in  reducing  the  flow  of  drugs  into 
this  country  is  only  one  part  of  the  overall  struggle,  but,  as  Secretary  Cheney  said,  the 
military  can  have  a  great  impact  in  this  effort. 


B.       WHAT  IS  DOD'S  INVOLVEMENT  NOW? 


Since  the  Department  of  Defense's  first,  somewhat  reluctant,  involvement  in 
counter-narcotics  in  the  1970's,  the  military  has  been  called  upon  for  increasing  support 
each  year.  The  Defense  Authorization  Act  of  1989  presented  DoD  with  a  primary  role 
in  the  drug  war.  We  now  see  all  segments  of  the  military  accepting  counter-narcotics 
assignments,  including  operational  missions,  such  as  ship  and  aircraft  patrols  with 
embarked  law  enforcement  personnel,  intelligence  gathering  and  surveillance,  and  material 
support.  The  overall  DoD  chain  of  command  for  counter-narcotics  is  shown 
in  Figure  1. 
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Figure  1    DoD  Counter-narcotics  Chain  of  Command. 


Obviously  there  is  high-level  attention  being  paid  to  the  problem.  The  three  Joint 
Task  Forces  command  the  military  assets  involved  in  actual  day-to-day  interdiction 
operations. 

1.  Joint  Task  Force  4 

Located  in  Key  West,  Florida,  JTF4  is  in  charge  of  operations  in  the  Atlantic 
area  of  responsibility  (AOR).  As  the  first  Joint  Task  Force,  their  mission  included  the 
creation  of  a  joint  fusion  center  for  tactical  intelligence  and  the  communications  necessary 
to  assemble  the  data  and  provide  a  finished  product  to  all  users.  JTF4  conducts  their  own 
detection  and  monitoring  operations  in  the  Atlantic  AOR,  and  coordinates  the  operations 
of  other  agencies,  such  as  the  Coast  Guard,  Customs  Service  and  local  law  enforcement 
activities. 

Although  JTF4  has  no  dedicated  assets,  they  do  have  tactical  control  of 
various  personnel,  ships,  patrol  and  intercept  aircraft,  and  intelligence  collecting  assets 
assigned  from  other  DoD  commands,  the  Coast  Guard  and  Customs  Service. 
[Ref.  10:  p.  78] 

2.  Joint  Task  Force  5 

Located  in  Alameda,  California,  JTF5  has  a  similar  responsibility  and  mission 
in  the  Pacific  AOR  as  that  of  JTF4  in  the  Atlantic.  They  conduct  organic  operations  to 
detect  and  monitor  aircraft  and  surface  vessels  suspected  of  smuggling,  as  well  as 
coordinating  other  agencies'  operations  within  their  AOR.  JTF5's  mission  includes 
integrating  into  the  intelligence  communications  network  and  providing  drug-related 
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intelligence  to  law  enforcement  agencies  as  appropriate.  [Ref.  11]  As  is  the  case 
with  JTF4,  JTF5  has  no  dedicated  assets,  but  takes  tactical  control  of  various  assets 
assigned  from  DoD  commands,  the  Coast  Guard  and  Customs  Service. 

3.       Joint  Task  Force  6 

Located  at  Ft.  Bliss,  Texas,  JTF6  has  a  much  different  mission  from  that  of 
JTF's  4  and  5.  Specifically,  "Joint  Task  Force  Six  plans  and  coordinates  all  DoD  support 
requested  by  Federal,  State  and  local  law  enforcement  agencies  within  the  Southwest 
Border  region."  [Ref.  12] 

Requests  for  support  usually  come  through  Operation  Alliance,  an  agency 
comprised  of  representatives  from  Federal,  State,  and  local  law  enforcement  agencies  of 
the  Southwest  Border  states:  Texas,  New  Mexico,  Arizona,  and  California.  Support 
provided  includes  observation  posts,  reconnaissance,  terrain  denial  operations  designed 
to  deter  smugglers  from  using  areas  along  the  border,  training  teams,  military  drug 
detection  dogs,  and  construction.  Personnel  from  all  military  services  and  various 
National  Guard  units  have  been  called  upon  to  perform  these  functions. 

C.       WHERE  MAY  DOD'S  INVOLVEMENT  LEAD? 

There  are  as  many  different  opinions  concerning  DoD's  future  involvement  in  drug 
interdiction  as  there  are  people  contemplating  it.  These  opinions  range  from  complete 
disassociation  to  complete  involvement  in  all  aspects,  not  just  detection  and  monitoring. 
Many  rational,  and  some  not  so  rational,  arguments  are  given  for  both  extremes  and 
everything  in  between. 
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One  of  the  main  arguments  against  the  use  of  military  forces  in  counter-narcotics 
is  the  Posse  Comitatus  Act  (Title  18  U.S.  Code,  Section  1385),  which  prohibits  the  direct 
use  of  military  forces  for  civil  law  enforcement  activities  within  the  United  States. 
Because  of  this  law,  military  forces  involved  in  drug  interdiction  operations  may  not 
perform  search  and  seizure  or  make  arrests.  The  military  is  limited  to  its  monitoring  and 
detection  mission  and  support  of  law  enforcement  activities. 

Another  strong  argument  is  that  drug  interdiction  operations  have  a  detrimental 
effect  on  DoD's  primary  mission  to  maintain  a  military  capable  of  fighting  a  war  to 
defend  the  U.S.  or  its  interests.  This  argument  contends  that  the  fight  against  drugs  is  not 
truly  a  war,  and  that  performing  area  patrols  and  ship  boardings,  detection  and  intercept 
of  slow  moving  aircraft,  and  patrolling  a  border  area  are  not  suitable  training  operations 
for  the  military.  Following  the  end  of  the  Cold  War  (i.e.,  the  reduction  of  the  Soviet 
threat)  drug  smuggling  takes  its  place  among  "principle  threats  ...  that  replace  the  old 
Soviet      threat."      [Ref.  13:  p.  18]  This      means      that      drug      interdiction 

operations  need  no  longer  be  construed  merely  as  training  evolutions,  but  are  actually  part 
of  the  military's  primary  mission. 

A  third  argument  is  the  cost  of  military  drug  interdiction  operations.  It  was 
estimated  that  the  Air  Force  spent  approximately  $433,000  per  drug  bust  in  1986. 
[Ref.  14:  p.  7]  Many      believe      that     the      $1      billion-plus      military      drug 

interdiction  budget  could  be  better  spent  on  education  and  treatment  programs.  These 
programs  are  extremely  important,  since  interdiction  alone  cannot  halt  the  problem  of 
drug  abuse.     However,  education  and  treatment  are  more  effective  when  drugs  are  more 
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expensive  and  less  readily  available.  "Most  leaders  of  prevention  and  treatment  programs 
recognize  this;  their  task  is  made  easier  when  drug  enforcement  works." 
[Ref.  15:  p.  3]  The  goal  of  interdiction  is  to  "complement  and  support  our 
international  drug  control  activities  and  domestic  law  enforcement  programs,  ...  to  create 
an  integrated  supply  reduction  program."  [Ref.  15:  p.  65]  Interdiction  functions  not  only 
to  seize  drugs  destined  for  the  U.S.,  but  to  act  as  a  deterrent  to  potential  drug  smugglers, 
and  to  force  the  price  of  illegal  drugs  to  rise  through  increased  risks  the  smugglers  face. 
In  addition,  "Interdiction  provides  a  highly  visible  sign  to  other  nations  of  U.S.  interest 
in  reducing  drug  use."  [Ref.  3:  p.  2]  This  is  a  typical  extension  of  the  age-old  Navy 
mission  of  showing  the  flag.  Thus,  the  cost  per  seizure  of  military  interdiction  operations 
does  not  encompass  the  whole  military  contribution  to  the  overall  drug  control  strategy. 

Plainly,  the  military  has  an  important  role  to  play  in  combatting  the  drug  problem 
faced  by  the  United  States.  While  no  one  wants  to  see  the  military  take  on  civilian  law 
enforcement  roles,  DoD  has  the  manpower,  assets,  technology,  and  money  to  assist  and 
supplement  Federal,  State,  and  local  law  enforcement  agencies  in  reducing  the  supply  of 
drugs. 

While  DoD  maintains  its  role  in  detection  and  monitoring  of  illegal  drugs  entering 
the  country,  there  are  many  who  are  pushing  for  an  even  greater  military  involvement. 
In  1990,  Senator  Mitch  McConnel  (R-KY)  proposed  a  plan  for  shooting  down  suspected 
airborne  smugglers  who  refuse  to  land  for  inspection  [Ref.  16:  p.  A 19]. 
Although  this  plan  was  not  put  into  effect,  it  did  have  a  large  amount  of  support. 
Military  and  LEA  personnel  are  providing  training  and  technical  assistance  for  anti-drug 
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and  counterterrorist  forces  in  a  number  of  production  and  transit  countries.  Operation 
Bahamas,  the  Turks  and  Caicos  Islands  (OPBAT),  a  multinational  strike  force  in  place 
since  1983,  has  been  very  successful  in  preventing  smuggling  flights  into  the  Bahamas. 
It  is  certainly  conceivable  that  similar  strike  forces  could  be  established  with  other 
countries. 

The  completion  of  the  aerostat  balloon  network  and  the  Caribbean  Basin  Radar 
Network,  as  well  as  proposed  over-the-horizon  radar  systems,  would  significantly  increase 
the  ability  to  detect  and  monitor  air  traffic  in  smuggling  areas,  and  act  as  a  deterrent  to 
smugglers.  Aerostats  are  tethered  radar-equipped  balloons,  flown  at  about  10,000  feet, 
which  "provide  surveillance  of  oceanic  and  land  areas  right  down  to  the  surface." 
[Ref.  17:  p.  7],  with  a  range  of  approximately  200  nautical  miles.  The 
Caribbean  Basin  Radar  Network  is  a  series  of  ground  based  radars  designed  to  provide 
nearly  complete  coverage  of  the  Caribbean  basin.  Figure  2  displays  the  proposed 
coverage  of  the  aerostat  network  [Ref.  18:  p.  84]. 

RADM  Walter  Leland,  USCG,  said,  "We  know  that  when  the  aerostats  are  flying 
the  traffic  stops."  [Ref.  19:  p.  7]  However,  there  are  widely  differing 
reports  on  the  operational  reliability  of  the  aerostat  balloons.  Some  Air  Force  and 
Customs  Service  officials  call  the  system  too  expensive  and  fragile,  claiming  the  balloons 
are  not  operational  more  than  60%  of  the  time  because  of  severe  weather  and 
maintenance      problems      [Ref.  20:  p.  15].  At      the      same      time,      Pentagon 

figures  showed  the  balloons  to  be  operational  70-75%  of  the  time. 
[Ref.  21:  p.  S14]         Such     discrepancies     cast     doubt     on     the     effectiveness     of 
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Figure  2   Aerostat  Radar  Network. 


these  systems,  in  turn  making  it  harder  to  get  funding  to  continue  and  expand  the 
program. 

Finally,  there  have  been  a  number  of  popular  novels  recently  published  which 
provide  more  speculation  about  the  possibilities  of  the  military's  future  role  in  drug 
interdiction.  Most  notable  of  these  are  Tom  Clancy's  Clear  and  Present  Danger  and  Dale 
Brown's  Hammerheads. 

Although  many  widely  ranging  proposals  for  military  involvement  in  counter- 
narcotics  and  numerous  arguments  against  any  use  of  the  armed  forces  in  drug 
interdiction  exist,  with  the  current  trend  toward  downsizing  the  armed  forces  it  is  difficult, 
if  not  impossible,  to  predict  exactly  what  DoD's  future  role  in  the  nation's  drug  control 
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strategy  will  be.    It  is  fairly  certain  however,  that  the  U.S.  military  will  continue  to  be 

involved.    A  recent  report  by  Rep.  Les  Aspin  (D-Wis.),  Chairman  of  the  House  Armed 

Services  Committee,  foresees  the  continuation  of  a  military  role: 

The  drug  problem  must  be  attacked  from  both  sides  --  both  the  demand  for  drugs 
and  the  supply  of  drugs  must  be  reduced....  Although  the  United  States  will 
continue  to  rely  primarily  on  nonmilitary  means  in  its  effort  to  curtail  the 
international  trade  in  drugs,  Americans  will  also  want  U.S.  military  forces  engaged 
in  the  struggle  with  drug  traffickers,  primarily  in  an  interdiction  role. 
[Ref.  13:  p.  17] 
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ID.    SIMULATION  MODEL 

A.       THE  ORIGINAL  SOAR  MODEL. 

1.       Background. 

The  Simulation  of  Adaptive  Response  (SOAR)  model  is  part  of  a  RAND 
Corporation  study  on  the  drug  interdiction  program,  undertaken  at  the  request  of  the 
Under  Secretary  of  Defense  for  Policy  in  February  1986.  This  study  culminated  with  the 
publication  of  analysis  results  in  three  RAND  reports. 

The  primary  conclusion  of  RAND's  effort  was  that: 


Increasing  drug  interdiction  efforts  are  not  likely  to  greatly  affect  the  availability 
of  cocaine  in  the  United  States. ...This  conclusion  is  driven  primarily  by  the  small 
share  of  total  drug  distribution  costs  that  are  accounted  for  by  the  smuggling  sector. 
Only  about  10  percent  of  the  final  price  of  cocaine  comes  from  smuggling  costs  and 
profits.  [Ref.  3:  p.  xi] 


This  conclusion  does  not  take  into  account  other  reasons  for  military  involvement  in  drug 
interdiction,  such  as  those  discussed  in  the  previous  chapter.  It  does,  however,  suggest 
that  the  need  exists  for  a  tool  to  be  used  in  counter-narcotics  analysis.  RAND's  SOAR 
model  provides  a  reasonable  starting  point  for  developing  such  a  tool. 

2.       Model  Description. 

SOAR  is  a  dynamic  network  simulation  which  allows  drug  smugglers  a 
number  of  air,  sea,  and  land  smuggling  routes.  These  routes  are  generic  in  that  they  are 
not  associated  with  particular  geographic  routes.  The  route  used  for  a  particular  shipment 
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of  drugs  is  probabilistically  selected,  depending  on  the  costs  of  shipping  drugs  via  that 
route  and  the  expected  risk  costs  for  shipping  via  the  method  associated  with  that  route 
based  on  a  perceived  probability  of  the  shipment  being  seized.  This  probability  is 
estimated  by  a  time- weighted  history  of  unsuccessful  shipments  on  that  route,  and  will 
change  with  time  as  shipments  succeed  or  are  interdicted. 

Time  between  shipments  is  an  exponential  random  variable  with  an  input 
mean.  The  quantity  of  each  shipment  is  an  input  constant,  as  are  the  capacity  and  costs 
of  the  routes  and  shipping  methods.  SOAR  does  not  attempt  to  model  the  adaptation  of 
changing  shipment  sizes  in  response  to  interdiction  severity.  The  capacity  of  the  route 
and  the  size  of  the  shipment  are  used  to  determine  the  number  of  trips  that  each  shipment 
will  require.  These  trips  occur  on  the  same  route  at  the  same  time,  leading  to  a  saturation 
factor  which  increases  the  probability  that  the  shipment  is  seized. 

Interdiction  in  SOAR  is  modelled  by  a  general  probability  of  interdiction  for 
each  route.  "Military  resources  are  treated.. .simply  as  means  for  augmenting  particular 
activities,  thereby  raising  the  probability  of  success  in  those  activities."  [Ref.  3:  p.  8] 
This  probability  is  an  input  value  and  can  be  changed  for  predetermined  time  periods 
(phases).  When  multiple  trips  are  required,  the  probability  of  interdiction  on  the  route, 
for  that  shipment,  is  increased  as  a  result  of  the  saturation  factor. 

The  smugglers'  overall  goal  is  to  ship  as  large  a  quantity  of  drugs  as  possible 
at  their  perceived  lowest  cost.  Costs  to  the  smuggler  include  the  cost  of  the  drug, 
operating  costs  for  the  method  of  shipment,  the  cost  of  replacing  lost  assets  if  a  shipment 
is  seized,  and  pay  to  personnel.    Smuggling  personnel  are  assumed  to  be  risk  averse, 
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therefore  pay  is  assumed  to  change  as  the  square  of  the  perceived  risk  involved  in  using 
a  particular  route.  Increasing  interdiction  rates  will  affect  the  overall  cost  to  the 
smugglers  by  raising  each  of  these  costs. 

For  their  model,  the  RAND  authors  attempted  to  estimate: 

•  Quantities  shipped,  by  route,  in  a  given  year. 

•  Number  of  shipments,  by  route,  in  a  given  year. 

•  Number  of  vessels  identified  as  suspicious,  by  route,  in  a  given  year.  Of  those,  the 
number  pursued,  by  route  in  a  given  year.  Of  those,  the  number  resulting  in 
seizures,  by  route,  in  a  given  year. 

•  Estimates  of  the  compensation  resulting  from  the  likelihood  of  prison. 

•  Estimates  of  smugglers'  nonrisk  compensation  and  profits.  [Ref.  2:  p.  9] 

As  can  be  expected,  little  concrete  data  was  available,  resulting  in  the  use  of  informed 
guesses  for  many  of  the  input  parameters. 

Output  of  the  model  includes  an  echo  of  the  input  parameters  and  the  statistics 
described  in  Appendix  A,  Sections  A. 2.  and  B.l.  respectively.  The  output  is  an  average 
of  an  input  number  of  runs  desired.  Ten  runs  were  used  in  the  RAND  analysis,  providing 
"an  adequately  precise  estimate  of  the  overall  means,  especially  compared  with  the 
imprecision  of  some  of  the  input  data  and  some  of  the  assumptions  incorporated  in  the 
model."  [Ref.  2:  p.  24]  This  output  was  not  used  in  our  analysis,  but  is  included  in  the 
appendix  for  information. 

Computational  details  of  the  SOAR  model  are  described  in  Chapter  IV  of 
Simulation  of  Adaptive  Response,  and  are  reasonably  easy  to  interpret  from  the  model 
itself. 
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B.       THE  SOAR  SIMULATION  AS  THE  OBJECTIVE  FUNCTION  EVALUATOR. 

1.       Measures  of  Effectiveness. 

In  evaluating  drug  interdiction  efforts,  many  measures  of  effectiveness  have 
been  proposed.  Civilian  law  enforcement  agencies  are  concerned  with  numbers  of 
seizures  and  smugglers  prosecuted.  News  reports  regularly  include  the  street  value  of 
drugs  seized.  Because  of  its  unique  role  of  detection  and  monitoring  and  support  to  law 
enforcement,  and  the  fact  that  the  armed  forces  are  prohibited  from  direct  participation 
in  the  arrest  and  seizure  phases  of  interdiction,  the  military's  measure  of  effectiveness 
(MOE)  is  based  on  the  number  of  supporting  missions  provided.  This  MOE  is  not  very 
useful  for  a  model  which  attempts  to  find  optimal  locations  for  interdiction  assets. 
Although  it  indicates  the  level  of  military  support  to  reducing  the  flow  of  drugs  into  the 
United  States,  it  provides  no  quantitative  measure  of  the  effectiveness  of  the  overall  drug 
interdiction  program. 

The  quantity  of  drugs  seized  is  a  popular  MOE,  as  it  provides  an  evaluation 
of  how  much  drugs  are  prevented  from  entering  the  country.  However,  it  is  possible  that, 
when  interdiction  efforts  are  most  successful,  no  drugs  at  all  would  be  seized.  That  is, 
shipments  could  be  deterred,  smugglers  could  be  forced  to  use  the  most  expensive 
methods  of  shipping  to  avoid  interdiction,  etc. 

The  major  goal  of  drug  interdiction  is  to  reduce  the  consumption  of  illegal 
drugs.  Military  assets  are  an  integral  part  of  increasing  the  effectiveness  of  interdiction, 
thereby  helping  to  achieve  this  goal.  In  addition  to  seizing  drugs  and  equipment, 
interdiction  increases  the  risks  which  smugglers  face  and  causes  them  to  change  their 
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methods  of  operation.  All  of  these  effects  of  interdiction  result  in  higher  costs  to  the 
smugglers,  costs  which  are  eventually  passed  on  to  the  consumers  of  illegal  drugs. 
Although  some  users  are  likely  to  be  insensitive  to  increased  drug  prices,  "the  body  of 
economic  support  for  a  relationship  between  price  and  consumption  is  too  strong  to  allow 
much  doubt  that,  at  least  in  the  long  run,  higher  prices  will  lead  to  lower  consumption." 
[Ref.  3:  p.  20]  Therefore,  if  it  can  be  measured  or  estimated  with  some  confidence,  cost 
to  the  smugglers  is  an  MOE  which  reflects  the  overall  goal  of  the  drug  interdiction 
program. 

As  computed  by  the  SOAR  model,  total  cost  to  the  smugglers  includes  the 
costs  of  unsuccessful  shipments,  thereby  incorporating  the  MOE  of  the  amount  of  drugs 
seized.  Parameters  for  an  alternate  objective  function  value  could  easily  be  extracted 
from  the  simulation,  or  the  optimization  algorithm  could  be  applied  to  any  desired 
objective  function. 

The  simulation  used  as  the  objective  function  evaluator  for  the  optimization 
model  presented  in  this  thesis  is  a  slightly  modified  version  of  the  original  SOAR  model. 
The  changes  involve  the  attempt  to  model  smuggling  routes  with  a  geographical  structure 
and  calculating  the  probability  of  interdiction  based  on  the  location  of  military  assets. 

2.       Model  Assumptions. 

In  any  model,  many  simplifying  assumptions  are  required.  This  is  especially 
true  in  attempting  to  model  a  process  of  this  scale.  Perhaps  the  biggest  assumption  lies 
in  the  use  of  the  SOAR  model  outcome  as  the  objective  function  to  be  optimized.  That 
is,  we  assume  that  the  simulation  provides  a  reasonable  estimate  of  the  total  cost  to  the 
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smugglers.   This  includes  accepting  as  reasonable  the  assumptions  made  in  constructing 
SOAR,  which  are  described  in  the  RAND  publications. 

Another  major  assumption  involves  modelling  the  drug  smuggling  routes. 
Except  for  the  natural  chokepoints  in  the  Caribbean,  little  is  known  about  the  exact  routes 
being  used.  Therefore,  the  routes  presented  in  this  model  are  assumed  to  be 
representative  of  area  routes  of  unknown  width.  In  the  Mona  Passage,  Yucatan  Channel 
and  other  chokepoints,  the  width  of  sea  routes  is  constrained  by  the  presence  of  land  on 
either  side.  However,  in  the  Pacific,  Southwest  Border  area,  and  for  many  air  routes,  the 
width  could  be  hundreds  of  miles.  Assets  which  provide  coverage  of  a  route  are  assumed 
to  be  patrolling  the  area  which  that  route  represents.  Figure  3  displays  a  general  idea  of 
smuggling  routes  as  interpreted  by  the  Drug  Enforcement  Agency. 


«  *  New  Yort  City 


Figure  3   Drug  flows  from  South  America  into  the  United  States. 
Source:    DEA  Map  1  -  1987  (Revised  1989) 
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A  third  assumption  is  made  when  computing  the  probability  of  detection  on 
the  routes.  Each  asset  represents  a  sensor  which  adds  a  probability  of  interdiction  to  the 
route  that  it  covers.  It  is  also  assumed  that  an  asset  covering  a  leg  of  the  route  acts 
independently  of  other  assets  on  that  and  other  legs,  and  that  detection  on  any  leg  of  a 
route  can  be  taken  as  an  event  independent  of  all  other  legs  on  that  route.  Realistically, 
a  patrol  aircraft  on  a  route,  which  detects  a  suspicious  surface  contact  is  likely  to  be  in 
either  direct  or  indirect  contact  with  a  surface  asset  on  that  leg  or  a  subsequent  leg,  giving 
information  about  the  presence  of  that  contact  to  the  surface  asset  and  increasing  the 
surface  asset's  probability  of  detecting  and  interdicting  that  contact.  Similarly,  a  surface 
asset  may  detect  a  suspicious  aircraft  and  report  it,  etc.  For  modeling  purposes,  however, 
we  believe  these  assumptions  to  be  reasonable.  Independence  between  assets  and  routes 
will  result  in  SOAR  providing  an  upper  bound  on  the  measure  of  effectiveness.  This 
consequence  should  be  considered  when  estimating  the  probabilities  of  interdiction  input 
to  the  model  for  each  asset  type  against  each  method  of  smuggling. 

3.       Data  Structures  and  Computations. 

a.     Smuggling  Routes. 

The  representation  of  routes  is  accomplished  using  a  hierarchical 
adjacency  list  of  geographical  points.  The  geographic  points  are  listed  in  a  data  file 
similar  to  that  shown  in  Appendix  A  Section  A.5.,  which  includes  points  of  origin  for 
drug  shipments,  intermediate  route  points,  and  points  of  entry  into  the  United  States. 
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The  information  for  the  routes  themselves  is  taken  from  another  data  file. 
A  description  of  the  inputs  and  an  example  of  the  data  file  are  in  Appendix  A  Section 
A.3. 

b.  Assets. 

Appendix  A  Section  A. 4.  describes  the  input  information  and  a  data  file 
for  interdiction  assets.  Although  we  have  listed  assets  as  SHIP1,  PAIR1,  etc.,  asset  need 
not  necessarily  refer  to  a  single  ship,  aircraft  or  ground  unit.  An  asset  may  be  considered 
whatever  mix  of  elements  is  required  to  achieve  the  range  and  probability  of  interdiction 
input  to  the  model.  For  simplicity,  we  will  continue  to  refer  to  assets  as  single  units. 
Asset  types  are  used  to  describe  assets  with  the  same  range  and  probability  of  interdiction. 
By  assigning  additional  asset  types,  the  matrices  containing  this  information  could  be 
expanded  to  include  unique  input  data  for  particular  assets. 

c.  Computation  of  route  probabilities. 

A  shipment  of  drugs  will  be  interdicted  on  a  particular  route  with  a 
probability  computed  in  the  subroutine  PROBCOMP.  In  this  subroutine,  the  legs  of  each 
route  are  checked  to  determine  if  assets  are  located  such  that  the  leg  falls  within  the  range 
of  the  assets.  Up  to  this  range,  the  asset  contributes  a  constant  probability  of  interdiction 
to  the  route.  Beyond  this  range,  the  probability  decays  exponentially  as  a  function  of  the 
additional  distance.  The  probability  of  interdiction  on  a  route  is  then  computed  as  simply 
the  combination  of  independent  probabilities  for  all  assets  covering  all  legs  of  the  route. 
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The  probability  decay  function  chosen  results  in  probabilities  for 
individual  assets  decreasing  to  approximately  13.5%  of  the  original  value  at  twice  the 
asset  range.  For  example,  an  asset  with  a  given  constant  probability  of  0.10  to  100  miles 
will  have  a  probability  of  interdiction  equal  to  0.0135  at  200  miles.  This  method  of 
computing  probabilities  can  be  supported  by  considering  the  given  range  of  an  asset  as 
that  area  which  can  be  patrolled  with  a  constant  probability,  and  exponential  decay  of  the 
probability  representing  the  decreasing  effectiveness  of  the  patrolling  as  the  distance 
between  the  asset  and  the  route  increases  past  that  range. 
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IV.    THE  ANNEALING  MODEL 

A.       GENERAL  DESCRIPTION 

The  simulated  annealing  method  is  so-called  for  its  resemblance  to  the  annealing 
of  solids,  in  which  a  material  is  heated  to  the  point  where  all  particles  arrange  themselves 
randomly  in  the  liquid  phase,  and  is  slowly  cooled  so  that  all  particles  arrange  themselves 
in  the  low  energy  ground  state  of  a  corresponding  lattice.  [Ref.  22:  p.  2] 

The  method  starts  at  either  a  predetermined  or  random  initial  point.  In  the  feasible 
space  of  the  independent  variables,  a  random  walk  samples  the  objective  function 
[Ref.  23:  p.  210]  (analogous  to  the  energy  state  of  a  physical  system). 
Each  step  providing  an  improved  objective  function  value  is  automatically  accepted. 
Steps  which  are  detrimental  to  the  objective  function  value  are  accepted  according  to  a 
Boltzmann  probability  function,  p  =  exp(-AE/KBT),  which  is  dependent  on  the  energy  and 
temperature  of  the  physical  system.  This  allows  the  path  to  walk  out  of  local  and  global 
optimal  points  [Ref.  23:  p.  211].  The  process  is  often  called  a  Metropolis  loop,  named 
for  the  author  who  used  a  model  simulating  the  evolution  of  a  solid  to  thermal 
equilibrium  in  1953  [Ref.  24:  p.  14].  Obviously,  the  algorithm  could  run 
indefinitely,  walking  in  and  out  of  the  global  optima.  One  way  to  prevent  this  from 
happening  is  to  specify  the  number  of  steps  the  algorithm  is  to  take.  However,  this  in  no 
way  guarantees  reaching  the  global  optimum.  Recently  much  research  has  been  put  into 
developing  cooling  or  annealing  schedules.    Essentially,  such  schedules  decrease  the 
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probability  of  accepting  a  detrimental  step  by  decreasing  the  temperature  of  the  system 
as  the  global  optimum  point  is  approached.  The  random  walk  is  terminated  when  a 
predetermined  number  of  successive  trials  does  not  produce  an  acceptable  objective 
function  value.  [Ref.  23:  p.  212] 

B.  APPLICATIONS 

The  simulated  annealing  method  has  been  applied  to  optimization  problems  in 
numerous  fields,  including  physics,  operations  research,  numerical  analysis,  biology, 
materials  science,  game  theory,  code  design,  etc.  [Ref.  22:  p.  14]. 

There  are  many  advantages  and  disadvantages  associated  with  the  simulated 
annealing  method.  The  algorithm  is  applicable  to  many  different  optimization  problems 
and  is  easy  to  implement.  The  greatest  disadvantage  is  the  potentially  prohibitive  amount 
of  time  required  to  converge  to  a  near-optimal  solution.  This  is  dependent  on  the  cooling 
schedule  and  the  step  size  used  in  the  random  walk.  [Ref.  25:  p.  91] 

C.  THE  GENERALIZED  ANNEALING  ALGORITHM 

As  stated  previously,  simulated  annealing  algorithms  have  been  used  almost 
exclusively  in  applications  which  minimize  the  objective  function.  The  probability 
function  is  especially  well-suited  to  cases  where  the  objective  function  minimizes  to  zero, 
or  where  the  optimal  value  is  known.  In  these  cases  the  probability  of  accepting  a  non- 
improving  step  is  driven  to  zero  as  the  optimum  is  approached. 

In  the  case  where  the  optimal  value  is  not  known,  the  typical  algorithm  calls  for 
selecting  some  initial  estimate  of  the  optimal  value  ((j)m)  and  allow  the  random  walk  to 
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sample  the  objective  function  value  ((j)),  continuing  until  <J)  -  (j)  m  becomes  negative,  then 
decrease  (j)m  repeating  this  process  as  required.  On  the  other  hand,  if  the  initial  estimate 
of  (J)m  is  too  low,  it  is  increased  as  part  of  the  cooling  schedule,  so  that  (J)  -  <(>m  will 
approach  zero.  [Ref.  23:  p.  213] 

We  have  attempted  one  of  many  possible  applications  of  this  method  to  a 
maximization  problem,  by  taking  an  initial  guess  for  the  optimum  as  an  input  value,  as 
well  as  the  number  of  steps  to  be  taken  before  reducing  the  guess.  Cooling  the  guess 
(<|>m)  is  accomplished  by  keeping  track  of  the  highest  value  found  so  far  ((j>B),  and 
reducing  (J)m  according  to  the  function  <J)m+1  =  <|>m  -  a  *  (<J>m  -  <j>B),  where  a  is  a 
predetermined  fraction  which  controls  the  cooling  speed.  This  ensures  that  <j>m  is  never 
reduced  past  the  best  value  found  in  the  random  walk.  If  the  random  walk  produces  a 
value  which  is  greater  than  either  the  current  <j>m  or  the  input  guess,  then  both  are 
increased  a  significant  amount  above  the  higher  value  (we  have  chosen  <j)m  =  1.5  *  <j)B). 
The  cooling  of  (|)ni  again  proceeds  according  to  the  previously  described  method. 

In  our  application  of  the  algorithm,  we  begin  with  a  large  step  size  (Ar)  in  order  to 
sample  a  significant  portion  of  the  feasible  region.  As  §m  is  reduced  (or  cooled)  we  also 
reduce  the  step  size  to  a  minimum  value  (Armin)  by  a  similar  equation,  sampling  a  tighter 
area  of  the  region.  When  a  potential  point  of  convergence  is  found,  Ar  is  reduced  to  a 
small  value  (Arlow)  in  an  attempt  to  find  an  acceptable  move  in  the  very  close  vicinity. 
If  an  acceptable  move  is  found,  the  step  size  is  reset  to  the  previous  value  (Arold),  and  the 
annealing  continues. 
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This  algorithm  is  not  guaranteed  to  locate  the  objective  function  value  associated 
with  a  global  optimum.  The  value  (3  in  the  probability  function  controls  the  rejection 
ratio,  and  will  also  determine,  in  part,  how  close  to  the  optimal  value  the  algorithm 
converges.  For  very  low  values  of  (3,  the  algorithm  will  converge  at,  or  very  close  to,  the 
optimum,  when  the  optimum  value  is  known. 

The  feasible  space  and  step  generation  used  for  our  model  are  very  simple, 
consisting  of  upper  and  lower  limits  on  the  latitude  and  longitude  for  asset  locations,  and 
random  steps  with  size  Ar.  For  example,  we  did  not  attempt  to  prevent  assets  from 
moving  onto  land,  restricted  airspace,  or  within  search  range  of  other  assets.  These  are 
the  more  detailed  geographic,  operational,  and  political  constraints,  such  as  those 
discussed  in  Chapter  I,  Section  B,  which  could  be  used  to  control  possible  asset  moves. 
Such  constraints  could  be  represented  in  a  subroutine  that  tests  feasible  locations  and  a 
set  of  rules  controlling  the  possible  moves  of  assets  relative  to  each  other. 

Our  algorithm,  shown  in  Figure  4,  is  very  similar  to  that  presented  by  Bohachevsky, 
Johnson  and  Stein  for  a  minimization  problem  allocating  ballistic  missile  interceptors 
[Ref.  26].  The  following  terms  used  in  the  algorithm  are  defined: 

•  (j)  is  the  objective  function  value.  A  subscript  of  0  indicates  the  current  value,  1 
indicates  the  value  at  the  new  location.   Other  subscripts  are  as  defined  above. 

•  Q  is  the  feasible  space,  that  is,  the  set  of  latitudes  and  longitudes  where  assets  may 
be  located. 

•  x  is  the  set  of  latitudes  and  longitudes  for  all  mobile  assets.    A  subscript  of  0 
indicates  the  current  location,  1  indicates  the  new  location. 

•  Ar  is  the  step  size.  Subscripts  min  and  low  indicate  the  minimum  step  size  used 
when  cooling  and  the  step  size  desired  in  attempt  to  find  a  more  precise  solution. 
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Step  1:  Compute  the  objective  function  value  for  initial  locations,  <J>0. 
Set  <t>m  =  input  guess. 

Step  2:  For  each  mobile  asset,  generate  two  independent  standard  normal  random  variates, 
A,  and  A2,  and  compute  the  components  of  the  random  direction  vector  d: 
d,  =  A/(A/*2  +  A2**2)**0.5,  i  =  1,2. 

Step  3:  Compute  candidate  locations  for  assets,  set  x,  =  Xq  +  Ar  *  d. 

Step  4:  If  the  candidate  locations  are  infeasible,  ie,  x  €  Q,  repeat  Step  2. 

Step  5:  Compute  the  new  objective  function  value,  4>,,  and  set  A<j>  =  (J>1  -  <|>0. 

Step  6:  If  the  step  increases  the  objective  function  value  (A(J)  >  0),  accept: 
set  k  =  0,  Xg  =  x,,  <\>0  =  <\>v  and  n  =  n  +  1 . 
If  the  step  size  has  been  reduced  in  Step  8  (k  >  S),  reset  step  size  Ar  =  Arold. 

If  <t>1   >  <t>B'  Set  ^B  =  01- 

6a:      Increase  §m  and  initial  guess  if  new  objective  function  value  is  greater,  ie, 
If  <)>,  >  <})m,  set<t>m=  1.5  *<!>,. 
If  <)),  >  initial  guess,  set  guess  =  1 .5  *  <J>,. 
Set  n  =  0. 

6b:  If  N  accepted  steps  taken  (n  >  N),  cool  §m  and  the  step  size: 
set  0m  =  <j>m  -  a  *  (<|>m  -  <J>B)  and  Ar  =  Ar  -  a  *  (Ar  -  Armin). 
Set  n  =  0. 

Go  to  Step  2. 

Step  7:  If  the  step  is  a  decrease  (A(J)  <  0),  compute  probability  of  acceptance  and 
determine  whether  to  accept  or  reject. 

Set  p  =  exp(  p  *  A<(>  /(<t>m-<J>0))-   Generate  a  uniform  random  variate  V. 
If  V  >  p,  reject,  set  k  =  k  +  1.  If  V  <  p,  accept,  set  k  =  0,  Xq  =  X,,  <J>0  =  <J>r 

Step  8:  If  the  stopping  rule  has  not  been  reached  (k  <  S),  go  to  Step  2, 

otherwise,  if  k  =  S,  attempt  to  find  a  more  precise  solution  with  small  step  size, 
save  current  step  size,  Arold  =  Ar,  set  Ar  =  Arlow,  and  go  to  Step  2. 

Step  9:  If  the  final  stopping  rule  has  not  been  reached  (k  <  2  *  S),  go  to  Step  2, 
otherwise  stop  and  return  <))0  as  the  optimal  solution  found. 


Figure  4     Annealing  algorithm  to  optimize  interdiction  asset  locations  relative  to  a 
stochastic  network. 
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(3  is  the  control  parameter  for  the  probability  function  (analogous  to  Boltzmann's 
constant). 

a  is  the  fraction  for  reducing  tym  and  Ar 

N  is  the  number  of  accepted  steps  at  each  value  for  <j)m 

S  is  the  stopping  rule,  the  control  for  the  number  of  attempted  moves  rejected 
before  the  solution  is  declared  found. 

k  is  used  to  count  the  number  of  consecutive  rejected  moves 

n  is  used  to  count  the  number  of  accepted  moves  for  each  level  of  ()>m 

The  FORTRAN  code  for  the  model  can  be  found  in  Appendix  B. 


D.       RESULTS 

The  optimization  model  developed  in  this  thesis  was  applied  to  a  smuggling 
network  of  eleven  air  and  sea  routes  in  the  Caribbean  region.  Figure  5  displays  this 
network,  formulated  from  the  sample  data  sets  in  Appendix  A.  Some  areas  of  the  routes 
used  are  intuitive,  for  example,  those  running  through  the  Yucatan  Channel  and  the  Mona 
and  Windward  Passages.  However,  these  are  meant  to  be  illustrative  of  actual  routes,  and 
are  not  based  on  intelligence  or  other  official  information.  The  model  was  run  on  a 
number  of  different  computers,  SUN  SPARC  1+  and  2  workstations,  Tektronix  4336, 
Silicon  Graphics  IRIX,  and  an  AMDAHL  5990-500  mainframe  computer.  The  four 
UNIX  systems  yielded  identical  results,  while  the  VMS-based  mainframe  computer  found 
similar  results  to  the  UNIX  systems. 

Two  applications  of  the  model  were  tested.  The  first  examined  the  locations  at 
which  the  annealing  algorithm  converged.    Ten  runs  each  were  performed  with  five 
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Figure  5    Sample  Caribbean  smuggling  network 

mobile  assets,  three  ships  and  two  aircraft.  The  resulting  locations  from  one  run  arc 
plotted  in  Figure  6.  It  is  clear  that  some  of  the  locations  could  easily  be  improved.  This 
is  due  in  part  to  the  value  of  the  probability  function  control  parameter,  p,  used  to  ensure 
convergence  of  the  algorithm  in  a  reasonable  amount  of  time.  Each  of  the  ten  runs  found 
different  final  locations,  but  all  appeared  reasonable  and  the  optimal  values  found  were 
similar,  with  mean  equal  to  75153.2  (dollars  x  1000)  and  a  standard  deviation  of  7204.1. 
The  second  application  explored  the  response  of  the  model  to  a  changing  number 
of  assets.  Each  run  used  equal  numbers  of  ships  and  aircraft.  Parameters  of  assets  are 
from  the  sample  data  set  in  Appendix  A.  We  found  that  the  function, 
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Figure  6   Results  of  five  asset  run  (one  of  ten) 

(J)  =  615430.6  -  71.07  *  (80  -  ASSETS)  **  2.07,  determined  using  non-linear  regression, 
fit  the  data  well  (R2  =  0.9922).  This  particular  function  is  dependent  on  the  number  of 
assets  at  the  last  data  point,  and  is  extremely  simplistic,  with  total  numbers  of  assets  as 
the  independent  variable.  Runs  with  more  than  80  assets  will  change  the  parameters  of 
this  function,  but  will  not  significantly  increase  the  maximum.  Although  the  fit  of  the 
curve  is  good,  the  variance  at  each  level  of  assets  cannot  be  assumed  equal,  mainly  due 
to  the  small  sample  size  at  the  levels  (2). 

The  regression  model  itself  is  not  a  fundamental  result,  but  can  be  used  to  obtain 
useful  information.  With  80  assets  we  are  very  nearly  saturating  the  routes,  obtaining 
little  additional  benefit  from  adding  assets.  We  will  assume  for  the  following  discussion 
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that  80  assets  induce  the  highest  possible  total  cost  to  the  smugglers.  In  Figure  7  we  have 
plotted  a  curve  depicting  the  fraction  of  maximum  value  achieved  versus  assets.  This 
fraction  represents  the  regression  function  value  divided  by  615430.6,  the  function 
evaluated  at  80  assets.  From  this,  we  see  that  50%  of  the  maximum  attainable  effect  on 
the  smugglers  is  obtained  with  24  assets,  75%  with  36  assets.  In  other  words,  75%  of  the 
maximum  effectiveness  is  achieved  with  fewer  than  half  the  maximum  number  of  assets. 
It  is  highly  unlikely  that  such  quantities  of  assets  would  be  available  for  employment  at 
any  particular  time.  However,  these  results  are  entirely  dependent  on  the  specific 
parameters  used  and  our  simple  asset  mix. 
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V.    CONCLUSIONS  AND  RECOMMENDATIONS 

A.       CONCLUSIONS 

In  this  thesis  we  have  demonstrated  that  simulated  annealing  can  be  employed  to 
find  near-optimal  locations  for  interdiction  assets  relative  to  a  network  with  stochastic 
flows  that  are  determined  by  a  simulation.  Our  approach  also  demonstrated  an  application 
of  this  type  of  model  to  explore  the  relationship  between  the  capabilities  of  interdiction 
assets,  their  number,  and  the  resulting  effect  on  the  overall  effectiveness  of  the 
interdiction  effort.  Although  our  model  uses  a  simple  representation  of  interdiction  assets, 
combinations  of  assets,  and  the  geographic,  operational  and  political  constraints  which 
would  normally  be  used  to  restrict  feasible  asset  locations  in  the  drug  interdiction 
problem,  we  believe  that  the  results  demonstrate  this  methodology's  viability  for 
development  into  the  type  of  analytic  tool  described  in  Chapter  I,  Section  B. 

The  primary  obstacle  to  further  use  of  this  methodology  in  counter-narcotics 
operational  analysis  and  planning  is  the  questionable  reliability  and  limited  amount  of 
information  available  on  smuggling  activities.  The  credibility  of  results  from  this  type 
of  model,  including  the  validity  of  SOAR,  is  highly  dependent  on  the  accuracy  of 
information  used  in  constructing  the  route  network,  and  assigning  costs  to  smugglers  and 
parameters  to  assets.  The  suspect  reliability  of  smuggling  data  will  make  it  difficult  to 
assess  this  credibility,  short  of  actually  employing  operational  assets  and  observing  the 
resulting  effects  on  smuggling  and  drug  use.    Even  this  level  of  verification  would  not 
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necessarily  provide  an  accurate  measure  of  credibility  since  the  targets  of  drug  interdiction 
are  the  only  true  source  of  the  information  required,  and  they  are  not  likely  to  come 
forward  with  such  details. 

B.       AREAS  FOR  FURTHER  RESEARCH 

The  next  logical  step  in  the  development  of  an  analytic  tool  useful  to  the  counter- 
narcotics  program  is  to  apply  the  methodology  explored  in  this  thesis,  or  the  model 
presented  itself,  to  an  expanded  problem.  Such  an  expansion  should  include  all  primary 
cocaine  smuggling  regions  (Atlantic/Caribbean,  Pacific,  and  Southwest  Border), 
constraints  reflecting  the  geographic  and  political  considerations  for  asset  locations,  and 
operational  inputs,  such  as  budget  feasible  asset  combinations  and  capability  restrictions 
of  assets.  These  constraints  could  easily  be  included  in  a  subroutine  which  tests  proposed 
asset  moves  for  feasibility.  Analysis  of  the  results  from  an  expanded  problem  will 
provide  a  more  realistic  evaluation  of  the  utility  that  may  be  derived  from  this  type  of 
tool. 

Expanding  the  size  of  the  problem  produces  a  major  drawback  in  the  use  of  the 
model  developed  in  this  thesis,  extensive  computer  run  times.  With  a  network  of  only 
eleven  routes,  the  five  asset  runs  described  in  our  results  averaged  seven  hours  each  on 
the  SUN  SPARC  1+  workstation.  Adding  routes  to  the  problem  will  significantly 
increase  the  run  time  of  SOAR.  Furthermore,  our  testing  used  only  one  trial  of  SOAR 
to  obtain  each  new  objective  function  value,  while  future  testing  should  use  the  average 
of  multiple  trials  to  reduce  the  variability  of  the  simulation  results.  This  will  increase  run 
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times  even  more,  by  approximately  a  multiple  of  the  number  of  trials.  For  these  reasons, 
analysis  of  the  SOAR  model  to  optimize  performance,  or  translation  into  a  more  efficient 
programming  language,  and  the  use  of  high  speed  computers  would  be  beneficial  to  future 
efforts. 

A  final  consideration  involves  the  variety  of  measures  of  effectiveness  (MOEs) 
being  used  by  drug  interdiction  agencies,  as  discussed  in  Chapter  in,  Section  B.l. 
Development  of  a  well-defined,  quantifiable  MOE  for  the  entire  drug  interdiction  program 
would  be  very  valuable  to  future  analysis.  The  MOE  produced  by  the  SOAR  model  for 
our  optimization,  total  cost  to  the  smugglers,  relies  heavily  on  information  which  may  not 
be  particularly  accurate,  as  would  be  the  case  for  any  MOE  that  reflects  the  goals  of  the 
interdiction  program  as  a  whole.  Continued  intelligence  efforts  to  gather  accurate 
information  on  illegal  drug  smuggling  activities  and  improvements  in  the  timeliness  of 
such  data  will  help  in  defining  an  overall  MOE  and  assessing  the  credibility  of  results 
from  analytic  tools  used  in  counter-narcotics  planning,  and  will  make  the  drug  interdiction 
program  more  effective. 
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A.    MODEL  INPUT  AND  OUTPUT 


INPUT  DATA  SETS 


1.        ANNEAL  DATA 


a.     Variables 


Variable  name 


STEPSIZE 


STEPSIZEMIN 


STEPSIZELOW 


Description 

Step  size  for  the  annealing 
algorithm  in  miles.  Real. 

Minimum  step  size  for  the 
algorithm  in  miles.  Real. 

Step  size  in  fine  search 
for  optimum  value  in  miles. 
Real. 


BETA 


Control  parameter  for 
annealing.  Adjusted  to 
control  the 

acceptance/rejection  ratio. 
Real. 


ALPHA 


The  stopping  rule  for  the 
algorithm.  This  number  of 
consecutive  rejections  stops 
the  algorithm.  Integer. 

Parameter  used  to  decrement 
the  estimate  of  the  optimal 
value  and  step  size. 
Adjusted  to  control  the  speed 
of  cooling.  Typically  between 
0.01  and  0.2.  Real. 
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GUESS 


Inital  estimate  of  the 
optimal  objective  function 
value.  Real. 


BIGN 


The  number  of  steps  to  be 
taken  at  each  level  of 
cooling.  Integer. 


b.     Sample  data  set 


INPUT  PARAMETERS  FOR  THE  ANNEALING  ALGORITHM 


STEPSIZE  (MILES) 
MIN  STEP  SIZE 
SMALL  STEP  SIZE 
CONTROL  PARAMETER 
STOPPING  RULE 
ALPHA 

GUESS  AT  MAXIMUM 
STEPS  PER  LEVEL 


7  5.0 
4  0.0 
5.0 
1.5 
100 
0.1 

100000000.00 
1000 


2.       SOAR  DATA 


a.     Variables 


Variable  name 
NUMMETHODS 

NUMPHASES 
NUMROUTES 

NUMTRIALS 


ENDTIME 


Description 

Number  of  methods  for 
smuggling  drugs.  Integer. 

Number  of  phases.  Integer. 

Number  of  routes  available. 
Integer. 

Number  of  times  the  analysis 
period  is  to  be  simulated. 
Reported  results  will  be 
averaged  over  the  trials. 
Integer. 

Number  of  days  to  be 
analyzed,  not  including  the 
runin  period. Integer. 
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RUNIN 

IX 

DRUGCOST 

EXSHIPMENTINTVL 

EXSHIPMENTSIZE 
METHODNAME(M) 

CAPCOST(M) 


CAPACITY(M) 


Number  of  days  to  run  the 
simulation  before  commeneing 
the  analysis.  Integer. 

Seed  for  random  number 
generator.  Integer. 

The  cost  at  the  source  of  a 
kilogram  of  cocaine.  Real. 

Expected  time  between 
shipments.  Real. 

Expected  shipment  size.  Real. 

The  name  of  drug  smuggling 
method  M.  Character. 

The  cost  to  the  smuggler  of 
an  unsuccessful  shipment 
using  method  M.  Does  not 
include  the  cost  of  the  drug 
or  the  cost  associated  with 
the  route.  Real. 

Maximum  amount  of  drug  that 
may  be  shipped  by  method  M. 
Real. 


b.     Sample  data  set 


METHODS 

PHASES    ROUTES    TRIALS   RUNIN 

ENDTIME 

MEMORY 

SEED 

3 

2         3 

2       120 

60 

0.1 

7243 

DRUGNAME 

DRUGCOST   SHI PMEAN 

SHIPSIZE 

DAILY  AMT 

COCAINE 

7500.0      0.71 

250.0 

350.0 

METHOD 

RISK  COMP  RISK  EXP 

SEIZE  COST 

MAX  SHIP 

AIR 

1200000.0       2.0 

200000.0 

2000  .0 

SEA 

1600000.0       2.0 

40000.0 

16000.0 

LAND 

10000.0       2.0 

5000.0 

50.0 

01 

61 
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3.        ROUTES  DATA 


a.     Variables. 


Variable  name 


ROUTENAME(R) 


ROUTECOST(R) 


ROUTEMETHOD(R) 


ROUTELOC(R) 


FIRST(R) 
NEXT(R,P) 


Description 

The  name  of  route  R. 
Character. 

The  cost  to  the  smuggler  for 
using  route  R,  incurred 
whether  or  not  the  shipment 
is  successful.  Real. 

The  index  of  the  smuggling 
method  that  is  used  on 
route  R.  Integer. 
(l=Air,2=Sea,3=Land). 

The  index  of  the  geographic 
area  of  route  R.  Integer. 
( 1  =Atlantic,2=Pacific, 
3=Southwest  Border). 

The  index  of  the  first  point 
of  route  R.  Integer. 

The  index  of  the  point 
following  P  on  route  R. 
Integer. 
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b.     Sample  data  set. 


ROUTENR   ROUTENAME  ROUTECOST    ROUTEMETHOD   ROUTELOC 
POINTS 

1  MOIMA  PASS  16000.00       2     1 

1  10  11    38 

2  WINDWARD  16000.00       2     1 

2  12  13    14    37 

3  YUCFLA1  16000.00       2     1 

3  15  36 

4  YUCALA1  16000.00       2     1 

3    15  35 


10 


11 


YUCTEX1 
3    15 

16000 

32 

.00 

YUCFLA2 
4    15 

160  0  0 
36 

.00 

YUCFLA3 
4    15 

16000 
35 

.00 

YUCTEX2 
4    15 

16000 
32 

.00 

CUBA 
3    36 

20000 

.00 

YUCALA2 
7    34 

20000 

.00 

YUCTEX3 

7    31 

20000. 

,00 
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4.        ASSET  DATA 

a.     Variables. 


Variable  name 
NUMASSETS 

NUMASSETTYPES 

ASSETNAME(A) 

LATASSET(A) 

LONASSET(A) 

ASSETTYPE(A) 


ASSETLOC(A) 


ASSETRANGE(AT,M) 


IPD(AT,M) 


Description 

Total  number  of  assets, 
mobile  and  fixed.  Integer. 

Total  number  of  asset  types. 
Integer. 

The  name  of  asset  A. 
Character. 

The  latitude  of  asset  A. 
Real. 


The  longitude  of  asset  A. 
Real. 


Index  of  the  type  of  asset  A, 

for  example, 

1  =air,2=ship,3=land, 


4=aerostat  balloon. 


Integer. 


Index  of  the  geographic- 
location  of  asset  A.  Same  as 
ROUTELOC.  Integer. 

Range  for  asset  type  AT 
versus  smuggling  method  M  for 
which  the  probability  of 
interdiction  is  a  constant. 
Real. 

Probability  of  interdiction 
for  asset  type  AT  versus 
smuggling  method  M  on  a 
route,  which  is  constant  out 
to  range  ASSETRANGE.  Real. 
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b.     Sample  data  set. 


TOTAL  NUMBER  OF  ASSETS  = 


NUMBER  OF  ASSET  TYPES  = 


PATROL  AIR, SHIPS, LAND  UNITS,  AEROSTATS 


lSSET 

NAME        ] 

LATITUDE 

LONGITUDE 

TYPE  ] 

jOC 

1 

PAIR1 

16.00 

75.85 

1 

1 

2 

PAIR2 

17.45 

80.95 

1 

1 

3 

SHIP1 

19.00 

68.00 

2 

1 

4 

SHIP2 

19.00 

75.00 

2 

1 

5 

SHIP3 

21  .98 

86.18 

2 

1 

6 

CUDLOI 

KEY  FL 

24.41 

81.30 

4 

1 

7 

MARFA  TX 

30.19 

104.01 

4 

3 

8 

GREAT  EXUMA 

23  .33 

75.47 

4 

1 

ASSET 

RANGE  1 

JS   METHOD 

METHOD 

ASSET 

AIR 

SEA 

LAND 

AIR 

100.0 

10  0.  0 

25.0 

SHIPS 

50.0 

50  .0 

0.0 

LAND 

10.0 

0.0 

5.0 

BALLOONS 

200.0 

0.0 

0.0 

P( INTERDICT) 

ASSET  VS 

METHOD 

METHOD 

ASSET 

AIR 

SEA 

LAND 

AIR 

0.  1 

0.1 

0.1 

SHIPS 

0.1 

0.1 

0.0 

LAND 

0.1 

0.0 

0.1 

BALLOONS 

0.1 

0.0 

0.0 

5.        GEOGRAPHIC  DATA 


a.     Variables. 


Variable  name 
POINTNAME(P) 

LAT(P) 

LON(P) 


Description 

The  name  of  geographic- 
point  P.  Character. 

The  latitude  of  point  P. 
Real. 

The  longitude  of  point  P. 
Real. 
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b.     Sample  data  set. 


POINT 

POINTNAME 

LATITUDE 

LONGITUDE 

1 

PT  GALLINAS 

12.15 

71.45 

2 

RIO  HACHA 

11.30 

7  3.00 

3 

SANTA  MARTA 

11  .15 

74  .10 

4 

TURBO 

8.00 

76.45 

5 

BUENAVENTURA 

3.58 

77  .50 

6 

TUMACO 

1.45 

78.45 

7 

MEDELLIN 

6.10 

75.47 

8 

BOGOTA 

4.32 

74.15 

9 

CALI 

3.20 

76.33 

10 

SMONAPASS 

17  .37 

68.00 

11 

NMONAPASS 

20  .50 

68.00 

12 

SWINDWARD 

17  .95 

75.18 

13 

NWINDWARD 

20.50 

73  .82 

14 

STRAIGHTS 

22.72 

76.91 

15 

YUCATAN 

21.98 

86.  18 

16 

CENT  MEXICO 

25.00 

103  .00 

17 

NOGOLES  AZ 

31  .20 

110.55 

18 

BISBEE  AZ 

31.30 

109.55 

19 

DOUGLAS  AZ 

31.20 

109.30 

20 

YSLETA  TX 

31  .42 

106.18 

21 

PRESIDIO  TX 

29.33 

104.23 

22 

DELRIO  TX 

29.21 

100.52 

23 

DOLORES  TX 

27  .42 

99.47 

24 

BROWNSVILLE  TX 

25.95 

97  .30 

25 

SAN  DIEGO  CA 

32  .41 

116.57 

26 

CALEXICO  CA 

32  .41 

115.30 

27 

??1        NM 

31.40 

108.50 

28 

??2       NM 

31.40 

108.25 

29 

FREEPORT  TX 

28.56 

95.21 

3  0 

GALVESTON  TX 

29.18 

94  .48 

31 

HOUSTON  TX 

29.46 

95.21 

32 

PORT  LAVACA  TX 

28.36 

96.38 

33 

SABINE  TX 

29.44 

93  .54 

34 

MOBILE  AL 

30.42 

88.03 

35 

PENSACOLA  FL 

30.25 

87  .13 

36 

ST  PETE  FL 

27.47 

82  .38 

37 

FLAMINGO  FL 

2  5.10 

80.55 

38 

DAYTONA  BCH  FL 

29.11 

81.02 

39 

MIAMI  FL 

25.45 

80.11 

4  0 

LOS  ANGELES  CA 

3  4.00 

118.15 

41 

BAJA  TIP 

22  .64 

109.82 

42 

S  MEXICO 

15.53 

96.18 

43 

PACIFICl 

13.32 

105.36 

44 

PACIFIC2 

30.00 

118.00 

45 

PACIFIC3 

20.40 

107  .27 

46 

PACIFIC4 

15.50 

105.00 
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B.       OUTPUT  DATA  SETS 

1.        SOAR  OUTPUT 

Output  of  the  model  includes  an  echo  of  the  input  parameters,  and  the 
following  statistics: 

1.  Expected  attempts  per  trial.  The  average  number  of  shipments  attempted  during  the 
period  being  analyzed. 

2.  Expected  successes  per  trial.  The  average  number  of  successful  shipments  during  the 
analysis  period. 

3.  Expected  interdictions  per  trial.  The  average  number  of  unsuccessful  shipments 
during  the  analysis  period.  Expected  attempts  =  Expected  successes  +  Expected 
interdictions. 

4.  Success  rate.   The  proportion  of  shipments  that  were  successful. 

5.  Interdiction  rate.   The  proportion  of  shipments  that  were  unsuccessful. 

6.  Cost  of  incomplete  shipments  (in  thousands).  The  average  cost  to  the  smuggler 
because  of  unsuccessful  shipments,  including  the  cost  of  the  method  (such  as  an 
airplane),  the  cost  of  the  drug,  and  the  cost  of  the  route  (such  as  gasoline). 

7.  Cost  of  completed  shipments  (in  thousands).  The  average  cost  to  the  smuggler 
because  of  successful  shipments,  including  the  cost  of  the  drug  and  the  cost  of  the 
route. 

8.  Total  cost  to  smugglers  (in  thousands).  The  sum  of  the  cost  of  incomplete  shipments 
and  the  cost  of  completed  shipments. 

9.  For  each  drug,  the  quantity  that  the  smugglers  attempted  to  ship. 

10.  For  each  drug,  the  quantity  that  the  smugglers  successfully  shipped. 

11.  For  each  drug,  the  quantity  that  the  smugglers  lost  because  of  unsuccessful  shipments. 

12.  For  each  route,  the  expected  attempts,  successes,  and  failures  are  reported. 

13.  For  each  phase  and  each  route,  the  expected  attempts,  successes  and  failures  are 
reported. 
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14.  For  each  drug  and  each  route,  the  expected  quantities  shipped,  captured,  etc.[Ref.  2: 
pp.  42,43] 

2.   ANNEALING  MODEL  OUTPUT 


MAX  TOTAL  COST  TO  SMUGGLERS  = 
(IN  THOUSANDS) 


78047 .95 


OPT  VALUE  FOUND  IN 


6630  STEPS 
831  STEPS  NOT  ACCEPTED 


419  5  INFEASIBLE  MOVES  FOUND 


FINAL  LOCATIONS 

ASSET  LATITUDE  LONGITUDE 


PAIR1 

22 

.715 

84  .834 

PAIR2 

23 

.539 

77 .752 

SHIP1 

20 

.728 

83.596 

SHIP2 

20 

.333 

86.613 

SHIP3 

21 

.687 

68.579 

CUDLOE  KEY  FL 

24 

.410 

81.300 

MARFA  TX 

3  0 

.190 

104.010 

GREAT  EXUMA 

23 

.330 

75.470 

FINAL 

PROBABILITY  OF 

INTERDICTION  BY 

ROUTE 

ROUTENAME 

PROBINTERDICT 

1 

MONA  PASS 

0.1774 

2 

WINDWARD 

0.2130 

3 

YUCFLA1 

0.3206 

4 

YUCALA1 

0.3176 

5 

YUCTEX1 

0.3176 

6 

YUCFLA2 

0.2910 

7 

YUCFLA3 

0.2879 

8 

YUCTEX2 

0.2879 

9 

CUBA 

0.2249 

10 

YUCALA2 

0  .2909 

11 

YUCTEX3 

0.1692 

CPU  TIME  =      17704.3  SECONDS 
MAX  VALUE  FOUND  =      7.83  917E+0  7 
FINAL  GUESS  VALUE  =      8.10845E+07 


NEXT  LOC  SEED  = 
RUN  COMPLETED 


298087429 
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APPENDIX  B.    FORTRAN  CODE 


PROGRAM  ANNEAL 

WRITTEN  BY  LT  JAMES  J.  HENRY  IV,  USN 
NAVAL  POSTGRADUATE  SCHOOL,  CODE  3  60 
OPERATIONS  ANALYSIS  CURRICULUM 


PROGRAM  VARIABLES 

INCLUDE  'COMMON. FOR' 

REAL  PROBMOVE, 
&      TLAT ( MAXASSETS )  , TLON ( MAXASSETS )  , 
&      OBJVAL0 , OBJVAL1 , OBJVALM, OBJVALB, DELOBJ , 
&      V, P,TOTCOST,CPUTIME, 
&      A(2) ,DENOM,TARRAY(2) 

INTEGER  K,LITTLEN, I , ROUTENR, MOVES, NOTMOVE, INFEAS 

LOGICAL  FEASIBLE, ACCEPT, FINISH, GOODINPUT 

CALL  SUBROUTINE  TO  READ  IN  THE  INITAL  DATA 


CALL  GETDATA (GOODINPUT) 

ISEED=IX 

ISEED1=IX 

IF (.NOT.  GOODINPUT)  STOP 

OBJVALM  =  GUESS 

DO  50  1=1, NUMASSETS 

TLAT ( I )  =  LATASSET ( I ) 
TLON ( I )  =  LONASSET ( I ) 
50  CONTINUE 

RUN  THE  SIMULATION  WITH  INITIAL  DATA  POINTS 

CALL  SMUGSIM (TLAT, TLON, TOTCOST) 

OBJVAL0  =  TOTCOST 

OBJVALB  =  OBJVAL0 

IF (OBJVALB  .GT.  GUESS)  GUESS  =  1.5  *  OBJVALB 

BEGIN  THE  ANNEALING 

MOVES  =  0 
NOTMOVE  =  0 
INFEAS  =  0 
FINISH  =  .FALSE. 
K  =  0 
LITTLEN  =  0 
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100  CONTINUE 

FIND  NEW  LOCATION  FOR  ALL  ASSETS  EXCEPT 
BALLOONS  (TYPE  4) 

DO  140  I  =  1,NUMASSETS 

IF(ASSETTYPEU)   .  NE .  4)  THEN 
120  CONTINUE 

CALL  LNORPCdSEEDl,  A,  2) 

DENOM  =  SQRT(A{1)**2  +  A(2)**2) 

TLAT(I)  =  LATASSET(I)  +  DELR*A ( 1 ) /DENOM 

TLON(I)  =  LONASSET(I)  +  DELR*A ( 2 ) /DENOM 

CHECK  FOR  FEASIBILITY 

IF(ASSETLOCd)  .  EQ .  1)  THEN 

CALL  LANTCHECK (TLAT ( I ) , TLON ( I ) , FEASIBLE] 
ELSEIF(ASSETLOC(I)   . EQ .  2)  THEN 

CALL  PACCHECK ( TLAT ( I ) , TLON {  I )  , FEAS I BLE ) 
ELSE 

CALL  SWBCHECK(TLATd)  ,  TLON  ( I )  ,  FEASIBLE) 
ENDIF 
IF d NOT.  FEASIBLE)  INFEAS  =  INFEAS  +  1 


IF d NOT.  FEASIBLE)  GO  TO  120 


14  0 


ENDIF 
CONTINUE 


RUN  THE  SIMULATION  WITH  NEW  DATA  POINTS 
CALL  SMUGS IM (TLAT, TLON, TOTCOST) 
OBJVAL1  =  TOTCOST 
DELOBJ  =  OBJVAL1  -  OBJVAL0 
ADJUST  MAX  VALUE  GUESS  IF  REQUIRED 
IF(OBJVALl  .GT.  OBJVALB)  OBJVALB  =  OBJVAL1 
IF(OBJVALB  .GT.  GUESS)  GUESS  =  1.5  *  OBJVALB 
IF (OBJVALB  .GT.  OBJVALM)  OBJVALM  =  GUESS 


TO  ACCEPT  THE  MOVE 


160 


DETERMINE  WHETHER  OR  NOT 

IF (DELOBJ  .GE.  0.0)  THEN 
ACCEPT  =  .TRUE. 

ELSE 

CK  =  OBJVALM  -  OBJVAL0 
IF(CK  .EQ.  0.0)  CK  =  0.000001 
P  =  EXP (BETA  *  DELOBJ  /  CK) 
CALL  LRNDPCdSEEDl,  V,  1  ) 
IF(P  .GE.  V)  THEN 


ACCEPT  =  .TRUE. 


ELSE 


ACCEPT  =  .FALSE. 

K  =  K  +  1 

NOTMOVE  =  NOTMOVE  +  1 
ENDIF 
ENDIF 

IF (ACCEPT)  THEN 
COUNT  =  0 
DO  160  I  =  1,NUMASSETS 

LATASSET(I)  =  TLAT ( I ) 

LONASSET(I)  =  TLON(I) 
CONTINUE 

MOVES  =  MOVES  +  1 
OBJVAL0  =  OBJVAL1 
LITTLEN  =  LITTLEN  +  1 
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ENDIF 

IF(K  .EQ.  S)THEN 

DELROLD  =  DELR 

DELR  =  DELRLOW 
ELSEIF(K  .GE.  (2  *  S) )  THEN 

FINISH  =  .TRUE. 
ELSEIF(K  .GT.  S  .AND.  ACCEPT)  THEN 

DELR  =  DELROLD 
ENDIF 

IF(LITTLEN  .GE.  BIGN)  THEN 

OBJVALM  =  OBJVALM  -  ALPHA  *  (OBJVALM  -  OBJVALB) 
DELR  =  DELR  -  ALPHA  *  {DELR  -  DELRMIN) 

ENDIF 

IF (OBJVALM  .LT.  OBJVALB)  OBJVALM  =  OBJVALB 


IF(.NOT.  FINISH)  GO  TO  100 
CLOSE (55) 

WRITE  THE  RESULTS 

OPEN(7  5,FILE='ANNEAL.OUT' ) 

WRITE(75, 520)  OBJVAL0/ 1000 . 0 

52  0  FORMAT (IX, 'MAX  TOTAL  COST  TO  SMUGGLERS  =  ',F10.2) 

WRITE(75, 530) 

53  0  FORMAT (1  OX,  '  (IN  THOUSANDS)') 

WRITE (75, 50  0) MOVES 

500  FORMAT (IX, 'OPT  VALUE  FOUND  IN  ',110,'  STEPS') 

WRITE ( 7  5 , 50 1 ) NOTMOVE 

501  FORMATQ9X,  110,  '  STEPS  NOT  ACCEPTED') 

WRITE(75, *) 
WRITE (75, 510) INFEAS 
510  FORMAT (IX, 110, '  INFEASIBLE  MOVES  FOUND') 
WRITE(75, *) 

WRITE(75, *) 

WRITE(75, *) 

WRITE(75, 534) 
534  FORMAT (1 IX, 'FINAL  LOCATIONS') 

WRITE(75, 535) 
53  5  FORMAT ( 1 IX, 'ASSET' , '  LATITUDE  ',' LONGITUDE  ') 

DO  550  I=1,NUMASSETS 

WRITE (75, 54  0) ASSETNAME ( I ) , LATASSET ( I ) , LONASSET ( I 
540     FORMAT (IX, A15.F10 .3.F10.3) 
550  CONTINUE 

WRITE(75, *) 

WRITE(75, 589) 
589  FORMAT( 'FINAL  PROBABILITY  OF  INTERDICTION  BY  ROUTE' 

WRITE (75, 59  0) 
59  0  FORMAT ( 'ROUTE'  , 5X,  ' ROUTENAME ' , 4X,  ' PROBINTERDICT '  ) 

DO  60  0  ROUTENR  =  1,  NUMROUTES 

WRITE (75, 610) ROUTENR, ROUTENAME (ROUTENR) , 
&  PROBINTERDICT (ROUTENR,  1) 

610     FORMAT (IX, 14, 5X, A10, 5X,F8.4) 
600  CONTINUE 
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CPUTIME    =    DTIME(TARRAY) 


WRITE(75, 

*  ) 

WRITE(75, 

*  ) 

WRITE(75, 

*  ) 

WRITE (75, 

*  ) 

WRITE (75, 

*  ) 

WRITE(75, 

*  ) 

WRITE(75, 

*  ) 

CLOSE (75) 

STOP 

END 

CPU    TIME    =     ', CPUTIME,'    SECONDS 
MAX    VALUE    FOUND    =     ' , MAXNOW 
FINAL   GUESS    VALUE    =     ' , OBJVALM 


NEXT    LOC    SEED 
RUN   COMPLETED' 


' , ISEED1 


******* 


******************* 


SUBROUTINES    TO    CHECK    FEASIBLE    REGIONS 

SUBROUTINE    LANTCHECK ( LAT , LON, FEASIBLE) 
* 

REAL    LAT , LON 

LOGICAL    FEASIBLE 
* 

FEASIBLE    =     .FALSE. 

IF (LON     .LE.     87.0     .AND.     LON     .GE.67.0     .AND. 

&  LAT     .LE.     25.0     .AND.     LAT     . GE .     14.0)     FEASIBLE    =     .TRUE. 

RETURN 

END 
********************************************************************* 


SUBROUTINE    PACCHECK (LAT, LON, FEASIBLE; 

REAL    LAT , LON 
LOGICAL    FEASIBLE 


FEASIBLE 


FALSE 


IF (LAT     .LE.     0. 

FEASIBLE    = 
ELSEIF(LON     . LE 

IF (LAT    .LE. 
ELSEIF(LON    . GE 

IF (LAT     .LE. 
ELSEIF(LON    .GE 

IF (LAT     .LE. 
ELSEIF(LON    .GE 

IF (LAT     .LE. 
ELSEIF(LON    .GE 

I F ( LAT 
ELSEIF(LON 

IF (LAT 
ELSEIF(LON 

IF (LAT 


LE. 
.GE 
LE. 
.GE 
LE. 
LAT 


0)    THEN 
. FALSE . 

.     120.0     .AND.    LON    .GE.     118.691)    THEN 
33.773)    FEASIBLE    =    .TRUE. 
.     118.0)    THEN 

32.50)     FEASIBLE    =     .TRUE. 
.     115.50)    THEN 

-132.70    -    1.40*LON))     FEASIBLE    =     .TRUE. 

113.00)    THEN 

-64.40    +     .80*LON))     FEASIBLE    =     .TRUE. 

107.00)     THEN 

-69.1667    +     .  8333*LON))    FEASIBLE    =     .TRUE 

84.50)    THEN 

-39.4444    +     .5556*LON))     FEASIBLE    =     .TRUE 

7  9.00)    THEN 

-15.5455    +     .2727*LON)        .AND. 

GE.     (121.5    -    1.5*LON))    FEASIBLE    =     .TRUE 


END  IF 

RETURN 

END 

•  it************************************************************- 


****** 


SUBROUTINE    SWBCHECK (LAT, LON, FEASIBLE) 

REAL    LAT , LON 
LOGICAL    FEASIBLE 
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FEASIBLE  =  .FALSE. 

* 

TF(LON  .LE.  117.13  .AND.  LON  .GE.  115.0}  THEN 

IF(LAT  .GE.  (33.74  -  ,0094*LON))  FEASIBLE  =  .TRUE. 
ELSEIF(LON  .GE.  111.0)  THEN 

IF(LAT  .GE.  (-7.015  +  .345*LON))  FEASIBLE  =  .TRUE. 
ELSEIFtLON  .GE.  108.295)  THEN 

IF(LAT  .GE.  (33.7421  -  .0222*LON))  FEASIBLE  =  .TRUE. 
ELSEIF(LON  .GE.  106.37)  THEN 

IF(LAT  .GE.  31.742)  FEASIBLE  =  .TRUE. 
ELSEIF(LON  .GE.  105.0)  THEN 

IF(LAT  .GE.  (-51.18  +  .7796*LON))  FEASIBLE  =  .TRUE. 
ELSEIF(LON  .GE.  104.474)  THEN 

IF(LAT  .GE.  (-173.537  +  1.9449*LON))  FEASIBLE  =  .TRUE. 
ELSEIF(LON  .GE.  103.474)  THEN 

IF(LAT  .GE.  (-38.3616  +  .6510*LON))  FEASIBLE  =  .TRUE. 
ELSEIFtLON  .GE.  103.36)  THEN 

IF(LAT  .GE.  29.0)  FEASIBLE  =  .TRUE. 
ELSEIFtLON  .GE.  102.731)  THEN 

IF(LAT  .GE.  (166.5411  -  1.3307*LON))  FEASIBLE  =  .TRUE. 
ELSEIF(LON  .GE.  101.462)  THEN 

IF(LAT  .GE.  29.837)  FEASIBLE  =  .TRUE. 
ELSEIF(LON  . GE .  100.747)  THEN 

IF(LAT  .GE.  (-72.6813  +  1.0098*LON))  FEASIBLE  =  .TRUE. 
ELSEIF(LON  . GE .  99.64)  THEN 

IF(LAT  .GE.  (-104.1222  +  1.3225*LON))  FEASIBLE  =  .TRUE. 
ELSEIF(LON  . GE .  99.00)  THEN 

IF(LAT  .GE.  (-168.6709  +  1.9703*LON))  FEASIBLE  =  .TRUE. 
ELSEIFtLON  .GE.  97.133)  THEN 

IFtLAT  .GE.  (1.7858  +  .2485*LON))  FEASIBLE  =  .TRUE. 
ENDIF 
RETURN 
END 

*******     SMUGSIM   SUBROUTINE     ************************************ 

SUBROUT INE  SMUGS IM ( TLAT , TLON , TOTCOST ) 
C        SMUGGLERS'  SIMULATION  -  2/3/87 

*  ADAPTED  FROM  RAND  CORPORATION  SIMULATION  OF  ADAPTIVE 

*  RESPONSE:  A  MODEL  OF  DRUG  INTERDICTION.  [Ref.  2] 

C  INSTALLATION  REMARKS: 

C 

C        YOU'LL  ALSO  NEED  TO  SELECT  "SIZING  PARAMETERS"  LARGE  ENOUGH 

C        TO  HANDLE  YOUR  ANALYSIS,  BUT  NOT  SO  LARGE  AS  TO  MAKE  THE 

C        PROGRAM  TOO  LARGE  TO  RUN  ON  YOUR  MACHINE.   SEE  THE  BEGINNING 

C        OF  COMMON. FOR  FOR  PARAMETERS. 

p*********************************************************************** 
C  PROGRAM  VARIABLES: 

INCLUDE  'COMMON. FOR' 

REAL  TLAT ( MAXASSETS ) , TLON ( MAXASSETS ) 

INTEGERM  T,  R 

p*********************************************************************** 
C  INITIALIZE    THE    DATA    COLLECTION    ARRAYS. 

DO    30,     R=l,MAXROUTES 
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AMOUNT ATTEMPTED (R)  =  0. 

AMOUNTSUCCEEDED(R)  =  0. 

ATTEMPTS (R)  =  0. 

SUCCESSES (R)  =  0. 

DO  20,  THISPHASE  =  1,MAXPHASES 

ATTEMPTSBYPHASE(R,THISPHASE)  =  0. 

SUCCESSESBYPHASE(R, THISPHASE)  =  0. 

2  0     CONTINUE 

3  0  CONTINUE 

SUCCESSCOSTS=0.0 
FAILURECOSTS=0 .0 

C         SET  UP  THE  TABLE  OF  WEIGHTS  FOR  PAST  SHIPMENTS. 
KMEMORY  =  ALOG(MEMORYVALUE) / (-LONGPAST) 
DO  40,  T=0,LONGPAST 

EXPTABLE(T)  =  EXP ( -KMEMORY*REAL (T ) ) 

4  0  CONTINUE 

*  CALL  THE  SUBROUTINE  TO  COMPUTE  ROUTE  PROBABILITIES  BASED  ON 

*  ASSET  LOCATIONS 

CALL  PROBCOMP(TLAT,TLON) 

C         CALL  THE  SUBROUTINE  THAT  DOES  THE  RUN  IN. 
CALL  PEACE 

C         DO  EACH  TRIAL. 

DO  50,  NTRIAL=1,NUMTRIALS 
CALL  SIMULATE 
50  CONTINUE 

*  COMPUTE  OUTPUT  TO  ANNEALING  ALGORITHM 

C        INITIALIZE  THE  TOTALS. 
TOTATTEMPTED=0. 
TOTSHIPPED=0. 
TOTATTEMPTS=0 . 
TOTSUCCESSES=0 . 
TOTSUCCESSCOST=0 . 
TOTFAILURECOST=0 . 

C        COMPUTE  THE  TOTAL  NUMBER  OF  ATTEMPTS,  SUCCESSES,  COSTS,  ETC 
DO  30,  R=l,NUMROUTES 

TOT ATTEMPTS  =  TOTATTEMPTS  +  ATTEMPTS (R) 

TOTSUCCESSES  =  TOTSUCCESSES  +  SUCCESSES (R) 

TOTATTEMPTED  =  TOTATTEMPTED  + 
&  AMOUNTATTEMPTED(R) 

TOTSHIPPED  =  TOTSHIPPED  +  AMOUNTSUCCEEDED ( R ) 

SUCCESSCOSTS  =  SUCCESSCOSTS  + 
&  AMOUNTSUCCEEDED (R)  *  DRUGCOST 

FAILURECOSTS  =  FAILURECOSTS  + 
&  ( AMOUNT ATTEMPTED { R ) - 

&  AMOUNTSUCCEEDED (R) )  * 

&  DRUGCOST 

3  0    CONTINUE 

TOTSUCCESSCOST=TOTSUCCESSCOST+SUCCESSCOSTS 

TOTFAILURECOST=TOTFAILURECOST+FAILURECOSTS 

C        GET  THE  AVERAGE  NUMBER  OF  ATTEMPTS,  SUCCESSES,  COSTS,  ETC. 
C        PER  TRIAL. 


53 


TOTATTEMPTS  =  TOTATTEMPTS/NUMTRIALS 
TOTSUCCESSES  =  TOTSUCCESSES/NUMTRIALS 
TOTSUCCESSCOST  =  TOTSUCCESSCOST/NUMTRIALS 
TOTFAILURECOST  =  TOTFAILURECOST/NUMTRIALS 
TOTCOSTS  =  TOTSUCCESSCOST  +  TOTFAILURECOST 

RETURN 
END 

Q*  ******************************************************************  * 
Q******************************************************************** 

SUBROUTINE  COMPFACT ( RNUMERATOR, RDENOMINATOR,  R) 
C        GET  THE  NUMERATOR  AND  DENOMINATOR  OF  THE  "R"  FACTOR,  WHICH 
C        WILL  BE  USED  TO  INCREASE  THE  PROBABILITY  OF  INTERDICTION 
C        ON  ROUTES  WITH  HIGHER  THAN  AVERAGE  TRAVEL. 

q*********************************************************************** 
C  PROGRAM  VARIABLES: 

INCLUDE  'COMMON. FOR' 

INTEGER*4  R, RNUMERATOR, RDENOMINATOR, S 


C 


*********************************************************************** 


RNUMERATOR  =  0 

RDENOMINATOR  =  0 

DO  10,  S=0,LONGPAST 

RDENOMINATOR  =  RDENOMINATOR  +  PASTSHIPMENTS ( R , S ) 
10    CONTINUE 

DO  20,  S=0, RECENT PAST 

RNUMERATOR  =  RNUMERATOR+ PASTSHIPMENTS ( R, S ) 
2  0    CONTINUE 

RETURN 

END 


C 


******************************************************************** 


SUBROUTINE  GETDATA (GOODINPUT ) 
C        READ  IN  THE  INPUT  DATA. 

C  PROGRAM  VARIABLES: 

INCLUDE  'COMMON. FOR' 

REAL  STEPSIZEM 

INTEGER*4    DAY1 , DAY2 , DAYN, M, N, TEMP, TEMP2  , 
&  P, ROUTENR, I, J, PTS(MAXPTS) , 

&  METHOD, 

&  POINT, TYPE, ASSET 

CHARACTER*15    PO I NTNAME (MAX POINTS) 

LOGICAL  GOODINPUT 

LOGICAL    TEMPIF 

p*********************************************************************** 

C        GOODINPUT  WILL  INDICATE  WHETHER  THE  DATA  WAS  CLEAN  AND  THE 
C        SIMULATION  SHOULD  BE  RUN. 
GOODINPUT  =  .TRUE. 
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OPEN ( 55 , FILE= ' ANNEAL . DAT ' ) 

* 

READ (55, *) 
READ (55, *) 

READ (55, 4010) STEPS I ZEM 
READ(55, 4010)STEPSIZEMIN 
READ(55, 4010)STEPSIZELOW 
4010  FORMAT (2 OX, F5 . 1 ) 
*      APPROXIMATION  OF  STEP  SIZE  IN  DEGREES 
DELR  =  STEPSIZEM/60.0 
DELRMIN  =  STEPSIZEMIN/60.0 

DELRLOW  =  STEPSIZELOW/60 .0 

* 

READ (55, 4 02  0) BETA 

4020  FORMAT (20X,F5. 2) 

* 

READ (55, 4030) S 

4030  FORMAT (20X, 15) 
* 

READ(55, 4040) ALPHA 

4040  FORMAT(20X,F5.3) 
* 

READ(55, 4050)  GUESS 

4050  FORMAT (20X,F15. 2) 

* 

READ (55, 4060)  BIGN 

4  060  FORMAT (2 OX,  110) 
* 

CLOSE (55) 


:        READ  IN  THE  OVERALL  SOAR  SIMULATION  DATA. 

OPEN (76,  FILE=' INPUT.DAT' ) 

READ(76, *) 

READ  (76,5010)  NUMMETHODS, NUMPHASES, NUMROUTES, 
&  NUMTRIALS,RUNIN,ENDTIME,MEMORYVALUE, IX 

5010  FORMAT  ( 6 ( 15 , 4X ) , F10 . 5 , 2X, 15 ) 
:        CHECK  WHETHER  TOO  MANY  METHODS  HAVE  BEEN  REQUESTED. 

TEMPIF  =  (NUMMETHODS  .GT.  MAXMETHODS) 

CALL  WRTERROR (TEMPIF, GOODINPUT) 

:        CHECK  WHETHER  TOO  MANY  PHASES  HAVE  BEEN  REQUESTED. 
TEMPIF  =  (NUMPHASES  . GT .  MAXPHASES) 
CALL  WRTERROR (TEMPIF, GOODINPUT) 

:        CHECK  WHETHER  TOO  MANY  ROUTES  HAVE  BEEN  REQUESTED. 
TEMPIF  =  (NUMROUTES  . GT .  MAXROUTES ) 
CALL  WRTERROR (TEMPIF, GOODINPUT) 

:        CHECK  WHETHER  TOO  MANY  DAYS  HAVE  BEEN  REQUESTED. 
TEMPIF  =  (ENDTIME  .GT.  MAXDAYS) 
CALL  WRTERROR (TEMPIF, GOODINPUT) 

:        MAKE  SURE  THE  INITIAL  SEED  IS  ODD.    (IN  CASE  WE  RUN  ON  A  SUN 
IX  =  IX  /  2 
IX  =  IX  *  2  +  1 

:        READ  IN  THE  DRUG  RELATED  DATA  RECORDS. 
READ (76, *) 
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READ  (76,5020)  DRUGNAME, DRUGCOST, EXSHIPMENTINTRVL, 
&  EXSHIPMENTSIZE,DAILYAMOUNT 

5020  FORMAT  (A10,4F10.5) 
C        READ  IN  THE  METHOD  RELATED  RECORDS. 
TEMP  =  NUMMETHODS 

IF  (TEMP  .GT.  MAXMETHODS)  TEMP  =  MAXMETHODS 
READ (76, *) 
DO  40,  M=1,TEMP 

READ  (76,5030)  METHODNAME (M)  , 
&  RISKCOMP (M) , RISKCOMPEXP (M) , CAPCOST (M) , 

&  CAPACITY (M) 

5030      FORMAT(A5, 5X,4F10.5) 
4  0  CONTINUE 

IF  (NUMMETHODS  .GT.  TEMP)  THEN 
DO  50,  M=TEMP+1, NUMMETHODS 
READ  (76,5030) 
50      CONTINUE 
ENDIF 

C        READ  IN  THE  DAYS  WHEN  EACH  PHASE  ENDS.   SET  UP  THE  VECTOR 
C        INDICATING  WHICH  PHASE  IS  IN  EFFECT  FOR  EACH  DAY. 
DAY1  =  0 

TEMP  =  NUMPHASES 

IF  (TEMP  .GT.  MAXPHASES)  TEMP  =  MAXPHASES 
DO  70,  N=1,TEMP 

READ  (7  6,5040)  DAY2 
5040        FORMAT  (15) 
TEMP2  =  DAY2 

IF  (TEMP2  .GT.  MAXDAYS)  TEMP2  =  MAXDAYS 
IF  (TEMP2  .GE.  DAY1)  THEN 

WRITE  (77,6016)  N,TEMP2 
6016  FORMAT  (8X, 'PHASE  ',12,'  LASTS  THROUGH  DAY  ',14) 

DO  60,  DAYN=DAY1,TEMP2 

CURRENTPHASE(DAYN)  =  N 

6  0  CONTINUE 

DAY1  =  TEMP2  +  1 
ENDIF 

7  0    CONTINUE 

CLOSE (76) 

*  READ  IN  THE  NETWORK  ROUTES  AND  ROUTE  DATA 

OPEN ( 9 6, F I LE=' ROUTES. DAT' ) 
READ (96, *) 
READ (96, *) 

DO  2  50  ROUTENR  =  1 , NUMROUTES 
DO  210  P  =  1,  MAXPOINTS 
NEXT (ROUTENR, P)  =  0 
210     CONTINUE 

READ (96,5075) ROUTENAME ( ROUTENR ) , ROUTECOST ( ROUTENR )  , 
&  ROUTEMETHOD ( ROUTENR ) , ROUTELOC ( ROUTENR ) 

507  5     FORMAT ( 9X, A10 , 3X, F10 .2, 2X, 15, 15) 

DO  22  0  I=1,MAXPTS 
PTS(I)=0 
220     CONTINUE 

READ(96, 5080) FIRST (ROUTENR) , ( PTS ( I ) , I=1,MAXPTS) 
5080     FORMAT(10X,9I5) 
P=FIRST (ROUTENR) 
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J=l 

IF(P  .NE. 

0) 

THEN 

NEXT { ROUT ENR, 

P)  = 

P  = 

NEXT(ROUTENR, 

P 

J  = 

J  + 

1 

GO  TO  24C 

) 

ENDIF 

CONTINUE 

CLOSE (96) 

240 

PTS  ( J 


250 


*  READ  IN  THE  ASSET  DATA 

OPEN (98, FILE=' ASSET.DAT' ) 
READ (98, 53  00)NUMASSETS 
READ (98, 5300)NUMASSETTYPES 
53  0  0  FORMAT ( 2 5X, 15) 
C        CHECK  WHETHER  TOO  MANY  ASSETS  HAVE  BEEN  REQUESTED. 
TEMPIF  =  (NUMASSETS  . GT .  MAXASSETS) 
CALL  WRTERROR (TEMPIF, GOODINPUT) 

C        CHECK  WHETHER  TOO  MANY  ASSET  TYPES  HAVE  BEEN  REQUESTED. 
TEMPIF  =  (NUMASSETTYPES  .GT.  MAXASSETTYPES ) 
CALL  WRTERROR (TEMPIF, GOODINPUT) 
READ (98, *) 
DO  3  50  ASSET  =  1,  NUMASSETS 

READ(98, 532  0) ASSETNAME( ASSET)  , 
t  LAT ASSET ( ASSET )  , LONASSET ( ASSET )  , 

&  ASSETTYPE( ASSET) , ASSETLOC (ASSET) 

53  2  0  FORMAT (5X,A1 5, 2 (F10 . 3 ) , 2 ( 15) } 
350  CONTINUE 
READ (98, *) 
READ (98, *) 
READ (98, *) 
READ (98, *) 
DO  360  TYPE=1, NUMASSETTYPES 

READ (98, 532  5) ( RANGE (TYPE, METHOD) , METHOD=l , NUMMETHODS ) 

53  2  5     FORMAT (15X, 3 (F5.1, 5X) ) 
3  60  CONTINUE 

READ (98, *) 
READ (98, *) 
READ (98, *) 
READ (98, *) 
DO  390  TYPE=1, 4 

READ (98, 533  0) ( I PD (TYPE, METHOD) ,METHOD=l , 3 ) 
5330     FORMAT(15X,3 (F5.3, 5X) ) 
390  CONTINUE 
CLOSE (98) 

*  READ  IN  THE  GEOGRAPHICAL  DATA 

OPEN (95, FILE=' GEOGR.DAT' ) 
READ (95, *) 
392  READ (95, 54  00 , END=3  95 ) POINT, POINTNAME ( POINT)  ,  LAT (POINT)  , LON ( POINT) 

54  0  0     FORMAT (IX, 15 , 5X, A15 , 5X, F7 . 3 , 5X, F7 . 3 ) 

GO  TO  392 
395  CONTINUE 
CLOSE (95) 

RETURN 

END 
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p******************************************************************** 
p******************************************************************** 

SUBROUTINE  INITSIM 
C         INITIALIZE  FOR  THE  CURRENT  SAMPLE  POINT. 

p*********************************************************************** 
C  PROGRAM  VARIABLES: 

INCLUDE  'COMMON. FOR' 
INTEGER*4  D, R 
REAL  RANDNM 

p*****************************************+***************************** 

C        SET  THE  SHIPMENT  COUNTER  TO  ZERO. 
NUMSHIPMENT  =  0 

C        DETERMINE  WHEN  THE  NEXT  SHIPMENT  WILL  OCCUR. 
CALL  LRNDPC (IX, RANDNM, 1) 
NEXTEVENT  =  ENDTIME  +1.0 

NEXTSHIPMENT  =  -  EXSHIPMENTINTRVL  *  ALOG ( RANDNM) 
IF  (NEXTSHIPMENT  . LT .  NEXTEVENT)  THEN 

NEXTEVENT  =  NEXTSHIPMENT 
ENDIF 

C        GET  THE  PAST  FROM  THE  LONG  RUN  IN. 
DO  30,  D=0,LONGPAST 

DO  20,  R=l,NUMROUTES 

PASTSHIPMENTS(R,D)  =  PEACESHIPMENTS (R, D) 
PASTFAILURES(R,D)  =  PEACEFAILURES (R, D) 

2  0        CONTINUE 

3  0    CONTINUE 

RETURN 

END 
p******************************************************************** 


C 


****************************************^ 


SUBROUTINE  PEACE 
C        DO  THE  INITIAL  RUN  IN. 

C  PROGRAM  VARIABLES: 

INCLUDE  ' COMMON . FOR ' 

REAL  AMOUNT, PROBCAUGHT, RFACTOR, TEMPTIME,  RANDNM,  ACTUALRISKCOMP 
INTEGER*4  D,R, RDENOMINATOR, RNUMERATOR, T, TRIPS 
LOGICAL  TRIPSUCCESS 

p*********************************************************************** 

C        DETERMINE  THE  TIME  OF  THE  FIRST  SHIPMENT. 
NEXTEVENT  =  RUNIN  +1.0 
CALL  LRNDPC ( IX, RANDNM, 1) 

NEXTSHIPMENT  =  -EXSHIPMENTINTRVL  *  ALOG (RANDNM) 
IF  (NEXTSHIPMENT  . LT .  NEXTEVENT)  THEN 

NEXTEVENT  =  NEXTSHIPMENT 
ENDIF 
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C        INITIALIZE  THE  ARRAYS  DESCRIBING  THE  PAST. 
DO  23  0,  D=0, LONG PAST 

DO  22  0,  R=l,NUMROUTES 

PASTSHIPMENTS(R, D)  =  0 
PASTFAILURES(R,D)  =  0 

22  0        CONTINUE 

23  0    CONTINUE 

C        FOR  EACH  DAY  OF  RUNIN. . . 
DO  100,  DAYNOW=l, RUNIN 
THISPHASE  =  1 

C  SHIFT  THE  ARRAYS  DESCRIBING  THE  PAST. 

DO  20,  R=l,NUMROUTES 

DO  10,  D=LONGPAST, 1, -1 

PASTSHIPMENTS(R,D)  =  PASTSHIPMENTS { R, D-l ) 
PASTFAILURES(R,D)  =  PASTFAILURES { R, D-l ) 
10  CONTINUE 

PASTSHIPMENTS (R, 0)  =  0 
PASTFAILURES (R, 0)  =  0 

2  0        CONTINUE 

C  LOOP  THROUGH  THE  DAYS  SHIPMENTS. 

3  0         IF  (NEXTEVENT  .GE.  DAYNOW  +  1.0)  GO  TO  100 

C  GET  THE  AMOUNT  OF  THE  NEXT  SHIPMENT. 

AMOUNT  =  EXSHIPMENTSIZE 

C  SELECT  THE  ROUTE  TO  BE  USED. 

CALL  SELROUTE ( R , TR I  PS , AMOUNT , ACTUALR I SKCOMP ) 

C  COMPUTE  THE  NUMERATOR  AND  DENOMINATOR  OF  THE 

C  "R"  FACTOR. 

CALL  COMPFACT ( RNUMERATOR , RDENOMINATOR , R ) 

C  FOR  EACH  TRIP  REQUIRED  TO  GET  AMOUNT  SHIPPED. . . 

DO  70,  T=l, TRIPS 

C  COMPUTE  THE  "R"  FACTOR. 

RFACTOR  =1.0 
IF  (RDENOMINATOR  .GT.  0.0)  THEN 

RFACTOR  =  RNUMERATOR/ RDENOMINATOR 
TEMPT IME  =  NEXTEVENT 

IF  (TEMPTIME  .GT.  RECENTPAST)  THEN 
IF  (TEMPTIME  .GT.  LONGPAST ) 
&  TEMPTIME  =  LONGPAST 

RFACTOR  =  RFACTOR  *  TEMPTIME/RECENTPAST 
ENDIF 
ENDIF 
IF  (RFACTOR  . LT .  1.0)  RFACTOR  =1.0 

C  COMPUTE  THE  PROBABILITY  OF  INTERDICTION. 

IF  (PROBINTERDICT(R, THISPHASE) .GE. .9999)  THEN 

PROBCAUGHT  =1.0 
ELSE 

PROBCAUGHT  =  1.0  -  EXP (RFACTOR *ALOG ( 1 . 0  - 
&  PROBINTERDICT (R, THISPHASE) ) ) 

ENDIF 

C  DETERMINE  WHETHER  THE  TRIP  WAS  SUCCESSFUL. 
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CALL  LRNDPC( IX, RANDNM, 1) 

TRIPSUCCESS  =  (RANDNM  .GE.  PROBCAUGHT ) 

C  DO  THE  BOOKKEEPING.   (NOT  VERY  EXTENSIVE 

C  DURING  THE  RUN  IN. ) 

PASTSHIPMENTS(R, 0)  =  PASTSHIPMENTS ( R, 0 )  +  1 
IF  (.NOT.  TRIPSUCCESS) 
&  PASTFAILURES (R, 0)  =  PASTFAILURES ( R, 0 )  +  1 

RNUMERATOR  =  RNUMERATOR  +  1 
RDENOMINATOR  =  RDENOMINATOR  +  1 
7  0  CONTINUE 

C  GET  THE  TIME  AND  TYPE  OF  THE  NEXT  SHIPMENT. 

CALL  LRNDPC( IX, RANDNM, 1) 
NEXTSHIPMENT  =  NEXTEVENT  - 
&  EXSHIPMENTINTRVL  *  ALOG ( RANDNM) 

NEXTEVENT  =  RUNIN  +1.0 
IF  (NEXTSHIPMENT  . LT .  NEXTEVENT)  THEN 

NEXTEVENT  =  NEXTSHIPMENT 
ENDIF 
GO  TO  3  0 
100    CONTINUE 

C        SAVE  THE  LAST  LONGPAST  DAYS  FOR  USE  INITIALIZING  EACH  TRIAL. 
DO  13  0,  D=0, LONGPAST 

DO  12  0,  R=l,NUMROUTES 

PEACESH I PMENTS ( R , D ) = PASTSH I PMENTS ( R ,  D ) 
PEACEFAILURES (R, D) = PASTFAILURES (R,  D) 

12  0        CONTINUE 

13  0    CONTINUE 

RETURN 
END 


SUBROUTINE  SELROUTE (RCHOSEN, TRI PS , AMOUNT, ACTUALRISKCOMP) 
C        SELECT  THE  ROUTE  FOR  THE  NEXT  SHIPMENT. 

C  PROGRAM  VARIABLES: 

INCLUDE  'COMMON. FOR' 

REAL  AMOUNT , CUMPROB , ROUTEPROB , COSTROUTE ( MAXROUTES }  , 
&        ACTUALRISKCOMP, PROBCAUGHT (MAXROUTES) , 

&        TEMPNUMERATOR (MAXROUTES ) , TEMPDENOMINATOR (MAXROUTES ) , 
&        TOTCOST, WEIGHT, TEMPRISKCOMP (MAXROUTES) 

INTEGER*4  RCHOSEN, TRIPS, RMETHOD, R, S 

p*********************************************************************** 

C        INITIALIZE  THE  NUMERATOR  AND  DENOMINATOR  FOR  EACH  ROUTE. 
DO  10,  R=l,NUMROUTES 

TEMPNUMERATOR (R)  =0.0 
TEMPDENOMINATOR (R)  =0.0 
COSTROUTE (R)  =0.0 
10  CONTINUE 

C        FOR  EACH  DAY  TO  BE  CONSIDERED,  ADD  IN  ITS  CONTRIBUTION. 
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DO  30,  S=0,LONGPAST 

WEIGHT  =  EXPTABLE(S) 
DO  20,  R=l,NUMROUTES 

TEMPDENOMINATOR ( R )  = 
&  TEMPDENOMINATOR (R)  +  WEIGHT*PASTSHIPMENTS ( R,  S ) 

TEMPNUMERATOR(R)  =  TEMPNUMERATOR (R )  + 
&  WEIGHT*PASTFAILURES{R,S) 

2  0     CONTINUE 

3  0  CONTINUE 

C        FOR  EACH  ROUTE,  COMPUTE  THE  PERCEIVED  PROBABILITY  OF  BEING 
C        CAPTURED  AND  HENCE  THE  EXPECTED  COST  OF  USING  THE  ROUTE. 
C        THE  PROBABILITY  A  ROUTE  WILL  BE  CHOSE  WILL  BE  PROPORTIONAL 
C        TO  THE  INVERSE  EXPECTED  COST  OF  USING  THE  ROUTE,  SO  GET  THE 
C        TOTAL  OF  THE  INVERSE  EXPECTED  COSTS  OF  USING  EACH  ROUTE. 
TOTCOST  =0.0 
DO  40,  R=l,NUMROUTES 

RMETHOD  =  ROUTEMETHOD(R) 

TRIPS  =  INT (AMOUNT /CAPACITY (RMETHOD)  +  0.999) 

IF  (TEMPDENOMINATOR (R)  . GT .  0.001)  THEN 

PROBCAUGHT ( R )  =  TEMPNUMERATOR ( R ) /TEMPDENOMINATOR ( R ) 
ELSE 

PROBCAUGHT (R)  =0.0 
END  IF 
IF  (RISKCOMPEXP( RMETHOD)  . LE .  0.0001)  THEN 

TEMPRISKCOMP ( R ) =RISKCOMP ( RMETHOD) 
ELSE 

TEMPRISKCOMP (R) =RISKCOMP (RMETHOD)  * 
&  ( (2*PROBCAUGHT(R) ) **RISKCOMPEXP ( RMETHOD) ) 

ENDIF 

COSTROUTE(R)  =  1 . 0 / (TRIPS* ( PROBCAUGHT (R ) *CAPCOST ( RMETHOD) + 
&  ROUTECOST(R)  +  TEMPRISKCOMP ( R )) + 

&  PROBCAUGHT ( R) *AMOUNT*DRUGCOST) 

TOTCOST  =  TOTCOST  +  COSTROUTE(R) 

4  0    CONTINUE 

C        NOW  CHOOSE  THE  ROUTE,  WHERE  THE  PROBABILITY  OF  CHOOSING  A 
C        GIVEN  ROUTE  IS  PROPORTIONAL  TO  THE  INVERSE  COST  OF  USING 
C        THAT  ROUTE. 

CALL  LRNDPC(IX,ROUTEPROB, 1) 
CUMPROB  =0.0 
R  =  0 
50  CONTINUE 

R  =  R  +  1 

CUMPROB  =  CUMPROB  +  COSTROUTE { R) /TOTCOST 
IF  (CUMPROB  .LT.  ROUTEPROB)  GO  TO  50 
RCHOSEN  =  R 

RMETHOD  =  ROUTEMETHOD( RCHOSEN) 

TRIPS  =  INT (AMOUNT/CAPACITY (RMETHOD)  +  0.999) 

ACTUALRISKCOMP=TEMPRISKCOMP (RCHOSEN) 

RETURN 
END 
p******************************************************************** 

p******************************************************************** 

SUBROUTINE  SIMULATE 
C        THIS  SUBROUTINE  PERFORMS  THE  SIMULATION  FOR  A  SINGLE  SAMPLE 
C        POINT. 
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p*  ***************  *  ******************  *  **********************************  * 

C  PROGRAM  VARIABLES: 

I NCLUDE  ' COMMON . FOR ' 

REAL  AMOUNT, PROBCAUGHT, RFACTOR, TEMPTIME, RANDNM, ACTUALRISKCOMP, 
£<        TIMESHIPPED(MAXSHIPMENTS) 

INTEGER* 4  D, R, RDENOMINATOR, RNUMERATOR  ,  T,  TRIPS, 
L  ROUTEUSED(MAXSHIPMENTS) 

LOGICAL  TRIPSUCCESS,SUCCESS(MAXSHIPMENTS) 

p*********************************************************************** 

C        INITIALIZE  FOR  THIS  TRIAL. 
CALL  INITSIM 

C        FOR  EACH  DAY . . . 

DO  100,  DAYNOW=0, ENDTIME 

C  GET  THE  POINTER  INTO  THE  ARRAY  OF  INTERDICTION  PROBABILITIES. 

THISPHASE  =  CURRENTPHASE(DAYNOW) 

C  SHIFT  THE  PAST  HISTORY  ARRAYS. 

DO  20,  R=l,NUMROUTES 

DO  10,  D=LONGPAST, 1, -1 

PASTSHIPMENTS(R,D)  =  PASTSHI PMENTS (R, D- 1 ) 
PASTFAILURES(R,D)  =  PASTFAILURES ( R , D- 1 } 
10  CONTINUE 

PASTSHIPMENTS(R, 0)  =  0 
PASTFAILURES (R, 0)  =  0 
2  0        CONTINUE 

C  LOOP  THROUGH  THE  EVENTS  THAT  HAPPEN  TODAY . 

30         IF  (NEXTEVENT  .GE.  DAYNOW  +  1.0)  GO  TO  100 

C  DETERMINE  THE  NEXT  TYPE  OF  DRUG  TO  BE  SHIPPED  AND 

C  HOW  MUCH . 

AMOUNT  =  EXSHIPMENTSIZE 

C  SELECT  THE  ROUTE. 

CALL  SELROUTE ( R , TR I  PS , AMOUNT , ACTU ALR I SKCOMP ) 

C  GET  THE  NUMERATOR  AND  DENOMINATOR  OF  THE  "R"  FACTOR. 

CALL  COMPFACT ( RNUMERATOR , RDENOMINATOR , R ) 

C  MULTIPLE  TRIPS  MAY  BE  REQUIRED  TO  SHIP  THE  GIVEN 

C  AMOUNT  OF  DRUG  ON  THE  SELECTED  ROUTE,  DUE  TO 

C  CAPACITY  CONSTRAINTS.   THE  VARIABLE  TRIPS  CONTAINS 

C  THE  NUMBER  OF  TRIPS  THAT  WILL  BE  REQUIRED.   I  EXPECT 

C  THAT  THIS  VARIABLE  WILL  USUALLY  BE  CONTAIN  A  ONE 

C  AND  THAT  AS  A  RESULT  DO-LOOP  7  0  WILL  USUALLY  BE 

C  EXECUTED  ONLY  ONCE. 
DO  70,  T=l, TRIPS 

C  COMPUTE  THE  "R"  FACTOR. 

RFACTOR  =1.0 
IF  (RDENOMINATOR  .GT.  0.0)  THEN 

RFACTOR  =  RNUMERATOR /RDENOMINATOR 
TEMPTIME  =  NEXTEVENT 

IF  (TEMPTIME  .GT.  RECENTPAST)  THEN 
IF  (TEMPTIME  . GT .  LONGPAST) 
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&  TEMPTIME  =  LONGPAST 

RFACTOR  =  RFACTOR  *  TEMPTIME/ RECENTPAST 
ENDIF 
END  IF 
IF  (RFACTOR  . LT .  1.0)  RFACTOR  =1.0 

C  COMPUTE  THE  PROBABILITY  OF  INTERDICTION. 

C  COMPUTE  THE  PROBABILITY  OF  INTERDICTION. 

IF  (PROBINTERDICT{R,THISPHASE) .GE. .9999)  THEN 

PROBCAUGHT  =1.0 
ELSE 

PROBCAUGHT  =  1.0  -  EXP (RFACTOR *ALOG ( 1 . 0  - 
&  PROBINTERDICT(R,THISPHASE) ) ) 

ENDIF 

C  DETERMINE  WHETHER  THE  TRIP  WAS  SUCCESSFUL. 

CALL  LRNDPC(IX,RANDNM, 1) 
TRIPSUCCESS  =  (RANDNM  . GE .  PROBCAUGHT) 

C  DO  THE  REQUIRED  BOOKKEEPING. 

NUMSHIPMENT  =  NUMSHIPMENT  +  1 
SUCCESS (NUMSH I PMENT)  =  TRIPSUCCESS 
TIMESHIPPED (NUMSHIPMENT)  =  NEXTEVENT 
ROUTEUSED( NUMSHIPMENT)  =  R 
ATTEMPTS (R)  =  ATTEMPTS (R)  +1.0 
AMOUNT ATTEMPTED ( R )  = 
&  AMOUNTATTEMPTED(R)  +  AMOUNT/TRIPS 

ATTEMPTSBYPHASE(R,THISPHASE)  = 
&  ATTEMPTSBYPHASE(R,THISPHASE) +1 

IF  (TRIPSUCCESS)  THEN 

SUCCESSCOSTS=SUCCESSCOSTS+ACTUALRISKCOMP+ 
&  ROUTECOST(R) 

SUCCESSES (R)  =  SUCCESSES (R)  +1.0 
AMOUNTSUCCEEDED(R)  = 
V  AMOUNTSUCCEEDED(R)  +  AMOUNT/TRIPS 

SUCCESSESBYPHASE(R,THISPHASE)  = 
&  SUCCESSESBYPHASE(R,THISPHASE)  +  1 

ELSE 

PASTFAILURES (R, 0)  =  PASTFAILURES ( R, 0 )  +  1 
FAILURECOSTS=FAILURECOSTS+ACTUALRISKCOMP+ 
&  ROUTECOST(R)  + 

&  CAPCOST ( ROUTEMETHOD ( R ) ) 

ENDIF 

PASTSHIPMENTS(R, 0)  =  PASTSHIPMENTS (R, 0 )  +  1 
RNUMERATOR  =  RNUMERATOR  +  1 
RDENOMINATOR  =  RDENOMINATOR  +  1 
7  0  CONTINUE 

C  DETERMINE  WHEN  THE  NEXT  SHIPMENT  OF  THIS  DRUG  WILL  BE 

CALL  LRNDPC ( IX , RANDNM , 1 ) 
NEXTSHIPMENT  =  NEXTEVENT  - 
&  EXSHIPMENTINTRVL  *  ALOG ( RANDNM) 

C  DETERMINE  WHAT  THE  NEXT  SHIPMENT  WILL  BE. 

NEXTEVENT  =  ENDTIME  +1.0 
IF  (NEXTSHIPMENT  . LT .  NEXTEVENT)  THEN 

NEXTEVENT  =  NEXTSHIPMENT 
ENDIF 
GO  TO  3  0 
100    CONTINUE 
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RETURN 

END 
p******************************************************************** 


C 


******************************************************************** 


SUBROUTINE  WRTERROR (ERRORCONDITION, GOODINPUT) 

C        THIS  SUBROUTINE  WRITES  THE  ERROR  MESSAGE,  IF  NEEDED. 

* 

p*********************************************************************** 
C  PROGRAM  VARIABLES: 

LOGICAL  ERRORCONDITION, GOODINPUT 
p**************************************************************r********* 

IF  (ERRORCONDITION)  THEN 

WRITE  (77,6001) 
6001        FORMAT  ('  * **ERROR*** ' ) 

GOODINPUT  =  .FALSE. 

ENDIF 

RETURN 

END 
(—******************************************************************** 


********************************************************************* 

SUBROUTINE  PROBCOMP (TLAT, TLON) 

*  THIS  SUBROUTINE  COMPUTES  PROBABILITIES  OF  INTERDICTION  ON 

*  ROUTES  AS  A  FUNCTION  OF  ASSET  DISTANCE  FROM  ROUTE. 

*  PROBABILITY  OF  INTERDICTION  IS  INPUT  TO  SMUGSIM. 
* 

.  INCLUDE  'COMMON. FOR' 

* 

REAL  Ml , M2 , Bl , B2 , TEMPPI , LATINT, LONINT,  DST,  DEND,  LENGTH, DIST, A, 
&      TLAT(MAXASSETS) , TLON (MAXASSETS) ,GSDIS,FCTR 
INTEGER  P,Q,R 


*  COMPUTE  PROBABILITY  OF  INTERDICTION  FOR  EACH  ROUTE 

DO  405  R=l,NUMROUTES 

DO  4  00  N=1,NUMPHASES 

PROBINTERDICT(R,N)  =0.0 
400     CONTINUE 
405  CONTINUE 

*  FOR  EACH  ROUTE  COMPUTE  EQUATIONS  FOR  ROUTE  SEGMENTS 

DO  440  R=l,NUMROUTES 
P=FIRST(R) 
Q=NEXT(R, P) 
410     CONTINUE 

IF(LON(P)  .NE.  LON(Q)   .AND. 
&  LAT(P)  .NE.  LAT(Q))THEN 

M1=(LAT(P)  -  LAT(Q) ) / 
&  ( LON ( P )  -  LON ( Q ) ) 

M2=-l .0/M1 

B1=LAT(P)  -  Ml*LON(P) 
ENDIF 
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TEMPPI=0 .0 

DO  42  0  ASSET=1,NUMASSETS 

IF(LON(P)  .NE.  LON(Q)    .AND. 
LAT(P)  .NE.  LAT(Q))THEN 

B2=LATASSET (ASSET)  -  M2*LONASSET 
LATINT=(B1*M2  -  B2*M1 ) / (M2-M1 ) 
LONINT=(Bl-B2) / (M2-M1) 
ELSEIF(LON(P)   . EQ .  LON(Q))  THEN 
LATINT=LATASSET ( ASSET ) 
LONINT=LON(P) 


ASSET ) 


ELSE 


LATINT=LAT(P) 
LONINT=LONASSET (ASSET) 


END  IF 


DETERMINE  CLOSEST  POINT  BETWEEN  ASSET  AND  ROUTE  SEGMENT 


DST=GSDIS(LON(P) , LAT(P) , 
LONINT,LATINT) 
DEND=GSDIS(LON(Q) ,LAT(Q) , 
LONINT,LATINT) 
LENGTH=GSDIS(LON(P) ,LAT(P) , 

LON(Q) , LAT(Q) ) 
IF(DST  .LE.  LENGTH  .AND.  DEND 
DIST=GSDIS (LONASSET (ASSET) 
LONINT,LATINT) 
ELSEIF(DST  . LT .  DEND)  THEN 

DIST=GSDIS (LONASSET (ASSET) , LATASSET (ASSET 
LON ( P ) , LAT ( P ) ) 
ELSE 

DIST=GSDIS ( LONASSET (ASSET 
LON ( Q ) , LAT ( Q ) ) 
END  IF 


LE.  LENGTH)  THEN 
LATASSET (ASSET) , 


LATASSET (ASSET) 


COMPUTE  PROBABILITIES  (INDEPENDENCE  ASSUMEDO 
IF(DIST  .LE. 
&  RANGE(ASSETTYPE(ASSET) , ROUTEMETHOD ( R ) ) )  THEN 

A  =  1.0  -  I PD(ASSETTYPE( ASSET ),  ROUTEMETHOD (R)  ) 

ELSE 

IF (RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD (R) )   .Eg.  0.0) 
&  RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD (R) )  =  0.0001 

FCTR=( RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD (R) ) -DIST) / 
L  RANGE (ASSETTYPE (ASSET) , ROUTEMETHOD ( R ) ) 

A  =  1.0  -  I PD( ASSETTYPE) ASSET) , ROUTEMETHOD (R) ) 
&  *  EXP(FCTR) 


42  0 


430 


440 


1.0  - 


A 


ENDIF 

TEMPPI 
CONTINUE 
DO  43  0  N=1,NUMPHASES 

PROBINTERDICT(R,N 


CONTINUE 

P  =  Q 

Q  =  NEXT(R, P) 
IF(Q  .NE.  0)  GO  TO 
CONTINUE 
RETURN 


1.0  -  TEMPPI) ) 


=  1.0  -(1.0  -  TEMPPI) * 
1.0  -  PROBINTERDICT(R, N) ) 


410 


END 


p************************************************************ 


******** 


REAL  FUNCTION  GSDIS (  XI,  Yl,  X2 ,  Y2 
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*  FUNCTION  COMPUTES  THE  GREAT  CIRCLE  DISTANCE  IN  MILES 

*  BETWEEN  TWO  POINTS  EXPRESSED  IN  LATITUDE  AND  LONGITUDE 

REALM        XI,Y1,X2,Y2 

REALM        RARC,RLAT1,RLAT2,  RLON1  ,  RLON2  ,  RPD 

* LOCAL  DATA  INITIALIZATION 

DATA  RPD  /0.1745329433017307E-01/ 

* GREAT  CIRCLE  DISTANCE 

* 

*  INPUTS  ARE  TWO  POINTS  EXPRESSED  AS  (W.  LONGITUDE,  N.  LATITUDE) 
* CONVERT  DEGREES  TO  RADIANS 

RLON1  =  RPD  *  XI 

RLAT1  =  RPD  *  Yl 

RLON2  =  RPD  *  X2 

RLAT2  =  RPD  *  Y2 

* SPHERICAL  ARC  LENGTH  IN  RADIANS 

ARGUMENT  =  (SIN ( RLAT1 ) *SIN ( RLAT2 ) 
&         +COS(RLATl) *COS(RLAT2) *COS ( RLON1 -RLON2 ) ) 

IF (ARGUMENT  . GT .  1.0)  ARGUMENT  =1.0 
IF ( ARGUMENT  . LT .  -1.0)  ARGUMENT  =  -1.0 
RARC  =  ACOSf ARGUMENT) 

* ARC  LENGTH  IN  STATUTE  MILES  WITH  ROUGH  CORRECTION  FOR  OBLATENESS 

*  OBLATENESS  CONTRIBUTES  LESS  THAN  0.2  5%  CORRECTION. 

*  IN  U.S.,  DARC  IS  APPROXIMATELY  3959.U*RARC 

* 

GSDIS  =  3959 . 0E0*RARC 

RETURN 
END 

*******     RANDOM  NUMBER  GENERATING  SUBROUTINES   ********************* 

C     ADAPTED  FROM  LEWIS, OR AV, AND  URIBE  ENHANCED  SIMULATION  AND 

C     STATISTICS  PACKAGE  [Ref.  27] 

C     THIS  PROGRAM  WILL  GENERATE  A  VECTOR  OF  NORMAL  RANDOM  VARIABLES 

C     ACCORDING  TO  THE  SINE-COSINE  METHOD 


C 


SUBROUTINE  LNORPC ( ISEED, A, N) 

INTEGER   N, I, IND, ISEED 

REAL   A(N) ,U(2) ,S, W,XSTAR 

DOUBLE  PRECISION   PI 

DATA  PI/3 .14159265358979D0/ 

IND  =  1 

DO  100  1=1, N 

IND  =  -IND 
IF  (IND.GE.0)  GOTO  2  0 
10         CALL  LRNDPC( ISEED, U, 2,  0,  0: 
S  =SQRT(-2*ALOG(U(l) ) ) 
W  =  2*PI*U(2) 
XSTAR  =  S*COS(W) 
A(I)  =  S*SIN(W) 
GOTO  100 
2  0         A(I)  =  XSTAR 
100      CONTINUE 
RETURN 
END 
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SUBROUTINE  LRNDPC  (ISEED,U,N) 


INTEGER 

REAL 

DOUBLE  PRECISION 

C      D31M1=2**31  -  1 

C      D31   =2**31 

DATA  D31M1/2147483647 .DO/ 
DATA  D31  /2147483648.D0/ 
DSEED  =  DBLE(ISEED) 


N,  I,ISEED 

U(N) 

D31M1,  DSEED,  D3 1 


DO  5  I  =  1 ,  N 

DSEED  =  DMOD(950706376.D0*DSEED,D31Ml) 
DSEED  =  DMOD ( 16807 . D0*DSEED, D3 1M1 ) 
U(I)  =  DSEED  /  D31 

CONTINUE 

ISEED  =  INT (DSEED) 

RETURN 

END 


******* 


COMMON . FOR 


SET  PAMETERS 
INTEGER* 4 
£< 
& 
& 


LONGPAST,  RECENTPAST,  MAXDAYS, 

MAXMETHODS,  MAXPHASES,  MAXROUTES,  MAXSHIPMENTS, 

MAXASSETS , MAXASSETTYPES , MAXPOINTS , 

MAXPTS 


PARAMETER 


;LONGPAST=120,  RECENTPAST=2  0, 
MAXDAYS=73  0, 

MAXMETHODS=3,  MAXPHASES=12  ,  MAXROUTES=100 , 
MAXSHIPMENTS=60  00,MAXASSETS=100/ 
MAXASSETTYPES=4 , MAXPOINTS=50 , MAXPTS=8 ) 


COMMON  VARIABLES 


REAL 
REAL 


VARIABLES 


AMOUNT ATTEMPTED ( MAXROUTES) , 

AMOUNTSUCCEEDED( MAXROUTES) , 

ATTEMPTS (MAXROUTES) , 

ATTEMPTSBY PHASE (MAXROUTES, MAXPHASES) 

CAPCOST (MAXMETHODS) , 

CAPACITY (MAXMETHODS) , 

DAILYAMOUNT, 

DRUGCOST , EX PT ABLE ( 0 : LONGPAST ) , 

EXSH I PMENTINTRVL , 

EXSHIPMENTSIZE, 

FAILURECOSTS , KMEMORY , MEMORYVALUE , 

NEXTEVENT , NEXTSH I PMENT , 

PROBINTERDICT (MAXROUTES, MAXMETHODS) , 

R I SKCOMP (MAXMETHODS) , 

RISKCOMPEXP (MAXMETHODS) , 

ROUTECOST (MAXROUTES) , 

SUCCESSCOSTS, 

SUCCESSES (MAXROUTES) , 
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&  SUCCESSESBY PHASE (MAXROUTES, MAXPHASES) , 

&  LAT( MAX POINTS) , LON (MAX POINTS ) , 

&  LAT ASSET { MAXASSETS ) , LONASSET ( MAXASSETS ) , 

&  DELR,DELRMIN,DELRLOW, 

&  BETA, ALPHA, RANGE (MAXASSETTYPES , MAXMETHODS ) , 

I  I PD ( MAXASSETT Y  PES , MAXMETHODS )  , GUESS 

INTEGER  VARIABLES: 

INTEGERM  CURRENT  PHASE  (  0  rMAXDAYS  )  ,  DAYNOW,  ENDTIME,  NTRIAL, 

&  NEXTEVENTTY  PE , NUMMETHODS , NUMPH ASES , 

&  NUMROUTES , NUMSH I PMENT , NUMTR I ALS , 

&  PASTSHIPMENTS(MAXROUTES, 0: LONG PAST) , 

&  PASTEAI LURES (MAXROUTES, 0:LONGPAST) , 

&  ROUTEMETHOD ( MAXROUTES ) ,RUNIN, 

&  THISPHASE, 

&  PEACESH I PMENTS (MAXROUTES, 0 : LONG PAST) , 

&  PEACEFAILURES (MAXROUTES, 0 :LONGPAST) , 

&  NUMASSETS , ROUTELOC ( MAXROUTES ) , 

&  ASSETLOC ( MAXASSETS ) , 

&  ASSETTYPE (MAXASSETS) , NUMASSETTYPES , 

&  S, NEXT (MAXROUTES, MAXPOINTS) , 

L  FIRST (MAXROUTES ) , BIGN 

INTEGER  VARIABLES: 

INTEGERM  IX,ISEED 

CHARACTER  VARIABLES: 

CHARACTER  *  1 0     METHODNAME ( MAXMETHODS ) 

CHARACTER* 12     ROUTENAME (MAXROUTES ) 
CHARACTER* 15     ASSETNAME (MAXASSETS ) 

COMMON  AMOUNTATTEMPTED, AMOUNTSUCCEEDED, ATTEMPTS, 

L  ATTEMPTSBY PHASE, 

&  CAPCOST, 

&  CAPACITY, 

&  DAILY AMOUNT, 

&  DRUGCOST , EXPTABLE , 

&  EXSHIPMENTINTRVL, 

&  EXSHIPMENTSIZE, 

&  FAILURECOSTS,KMEMORY,MEMORYVALUE, 

Sc  NEXTEVENT,NEXTSHIPMENT, 

&  PROBINTERDICT, 

&  RISKCOMP, 

&  RISKCOMPEXP, 

&  ROUTECOST, 

&  SUCCESSCOSTS, 

&  SUCCESSES, 

&  SUCCESSESBYPHASE, 

&  LAT,LON, 

&  LAT ASSET , LONASSET , 

&  DELR,DELRMIN, DELRLOW, 

Sc  BETA,  ALPHA,  RANGE,  IPD,  GUESS, 

&  CURRENTPHASE, DAYNOW, ENDTIME,  NTRIAL, 

&  NEXTEVENTTY PE, NUMMETHODS, NUMPHASES, 

&  NUMROUTES , NUMSH I PMENT , NUMTRI ALS , 

&  PASTSH I PMENTS, 

&  PASTFAILURES, 

&  ROUTEMETHOD, RUNIN, 

&  THISPHASE, 
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&  PEACESHIPMENTS, 

&  PEACEFAILURES, 

&  NUMASSETS , ROUTELOC , 

&  ASSETLOC, 

&  ASSETTY  PE , NUMASSETTY  PES , 

&  S, NEXT, FIRST, BIGN, 

&  IX.ISEED, 

&  METHODNAME, 

&  ROUTENAME, ASSETNAME 
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